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ABSTRACT 


This  thesis  deals  with  improving  the  miss  distance  of  a  missile,  with  imaging  seeker(s), 
by  utilizing  dynamic  image  processing.  In  an  encounter  with  a  missile,  a  target  tries  to 
avoid  the  missile  by  performing  an  evasive  maneuver  when  the  missile  is  at  a  relative  dis¬ 
tance  which  maximizes  the  miss  distance.  Dynamic  image  processing  permits  us  to  identify 
the  evasive  maneuver  of  the  target  by  estimating  its  acceleration  in  magnitude  and  direc¬ 
tion.  This  thesis  studies  methods  of  utilizing  this  additional  information  about  the  target’s 
behavior  in  order  to  improve  the  missile’s  performance.  First  the  proportional  navigation 
guidance  iaw  is  explored  in.  order  to  verify  its  advantages  and  weaknesses.  Then,  methods 
of  obtaining  the  time  dependent  3-D  movement  of  a  target  from  its  image  plane  feature 
point  correspondences  are  derived.  The  3-D  components  of  tire  target’s  acceleration  are  ob¬ 
tained  by  using  a  Kalman  filter.  Missiles  with  two  cameras,  one  camera  and  one  seeker  (ra¬ 
dar  or  IR),  and  only  one  camera  are  considered.  Methods  to  get  stereo  vision  by  using  the 
one  camera  plus  one  seeker  setup  and  the  single  camera  setup  are  proposed.  Advanced 
guidance  laws,  namely  advanced  proportional  navigation  and  optimal  guidance  are  de¬ 
rived,  for  a  3-D  environment.  A  three  dimensional  simulation  program  is  developed  using 


classical  proportional  navigation,  advanced  proportional  navigation,  and  optimal  guidance. 


Tlie  engagement  is  simulated  using  state  variable  design  and  the  performance  of  the  guid¬ 


ance  laws  is  compared. 
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I.  INTRODUCTION 


The  U.S,  as  a  result  of  the  highly  effective  kamikaze  attacks  during  World  War  n  on 
U.S  vessels,  initiated  the  development  of  the  first  tactical  missile  (Lark  guided  missile). 
Since  that  time,  proportional  navigation  guidance  has  been  used  in  virtually  all  the  world’s 
endoatmospheric  tactical  radar,  infrared(IR),  and  television(TV)  guided  missiles. 
Proportional  guidance  works  well  not  only  for  predictable  targets,  but  also  for  highly 
responsive  ones  (i.c.  targets  executing  evasive  maneuvers).  The  proportional  navigation 
guidance  technology  currently  in  use  appears  to  be  adequate,  if  the  effective  time  constant 
of  the  guidance  system  is  short  in  comparison  with  the  flight  time  and,  if  the  missile  has 
considerable  acceleration  advantage  over  the  target.  The  popularity  of  this  interceptor 
guidance  law  is  the  result  of  its  simplicity  of  implementation,  and  effectiveness.  Although 
proportional  navigation  was  apparently  known  by  the  Germans  during  World  War  n,  no 
applications  of  it  were  reported.  In  the  U.S,  this  guidance  law  was  studied  under  the' 
auspicious  of  the  U.S  Navy.  Proportional  navigation  was  originally  conceived  from 
physical  reasoning.  The  mathematical  derivation  of  the  “optimatility”  of  proportional 
navigation  came  more  than  20  years  later. 

This  research  develops  a  three  dimensional  missiie/target  simulation  using  three 
techniques  of  interceptor  guidance,  namely  classical  proportional  navigation,  augmented 
proportional  navigation  and  optimal  guidance.  The  primary  research  goal  is  to  improve  the 
miss  distance  of  a  missile  with  imaging  seeker(s)  by  utilizing  dynamic  image  processing. 
The  existent  dynamic  image  processing  algorithms  can  be  used  to  estimate  motion 
parameters  of  the  target  This  additional  information,  about  the  target  behavior,  will  be 
included  in  the  proportional  navigation  homing  loop  in  order  to  increase  the  missile 
percentage  of  kill  by  improving  the  final  miss  distance.  Information  about  the  target  motion 
is  especially  important  in  the  final  phase  of  the  engagement,  given  that  an  evasive 
maneuver  performed  by  the  target  creates  appreciable  miss  distance  that  may  preclude  a 
target  kill.  In  an  encounter  with  a  missile,  a  target  tries  to  avoid  the  missile  by  performing 
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a  evasive  maneuver  when  the  missile  is  at  a  relative  distance  that  maximizes  the  miss 
distance.  A  simulation  of  the  adjoint  model  of  the  linearized  homing  loop  permits  us  to 
obtain  miss  distance  projections  as  a  function  of  flight  time  or,  if  preferable,  as  a  function 
of  the  time  to  go.  The  target  can  induce  the  most  miss  distance  by  executing  an  evasive 
maneuver  at  a  short  time  to  go.  More  precisely,  the  optimal  evasion  from  the  target  “point 
of  view”  would  be  a  series  of  maneuvers  at  the  times  of  flight  that,  by  superposition, 
produce  the  most  miss  distance.  Estimating  the  target  maneuver  and  incorporating  this 
information  into  the  guidance  control  input  is  perhaps  the  difference  between  success  and 
failure. 

Chapter  n  introduces  the  idea  of  proportional  navigation  and  how  the  actual  guidance 
law  is  developed.  Chapter  HI  deals  with  estimating  the  target  motion  parameters  by  using 
two  perspective  views.  In  Chapter  IV,  the  augmented  proportional  navigation  and  optimal 
guidance  will  be  derived.  Also  in  this  chapter,  a  tridimensional  missile/target  simulation  is 
developed  using  classical  proportional  navigation.  Subsequently,  the  target’s  estimated 
motion  parameters  will  be  incorporated  in  the  tridimensional  engagement  by  using 
augmented  proportional  guidance  and  optimal  guidance.  Chapter  V  consists  of  actual 
simulation  results.  The  different  control  laws  will  be  tested  and  compared  by  producing 
miss  distance  projections  for  different  evasive  maneuvers.  Finally,  conclusions  and 
recommendations  follow  in  Chapter  VI.  All  computer  simulations  are  developed  using 
Matrix  Laboratory(MATLAB). 
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n.  FUNDAMENTALS  OF  TACTICAL  MISSILE  GUIDANCE 


A.  GENERAL 

Proportional  navigation  guidance  (PROPNAV)  commands  the  missile  to  turn  at  a  rate 
proportional  to  both  the  angular  velocity  of  the  line  of  sight  (LOS)  and  the  closing  velocity. 
The  constant  of  proportionality  is  a  unitless  designer  chosen  gain  (usually  in  the  range  3-5) 
known  as  the  effective  navigation  ratio/constant  Mathematically,  the  guidance  law  can  be 
stated  as 

um  =  NVc *■'  (Eq2.i) 

where  um  is  the  acceleration  command  which  is  petpendicular  to  the  instantaneous  LOS .  N 
is  the  effective  navigation  constant  Vc  is  the  closing  velocity  along  the  LOS,  and  A.  is  the 
LOS  angle  (in  rad).  The  overdot  indicates  the  time  derivative. 

If  the  navigation  ratio  is  greater  than  1,  the  missile  will  be  turning  faster  than  the  LOS, 
and  thus  the  missile  will  build  up  a  lead  angle  with  respect  to  the  line  of  sight.  For  a  constant 
velocity  missile  and  target  the  generation  of  this  lead  angle  can  put  the  missile  on  a  collision 
course  with  the  target  (zero  angular  velocity  of  the  line  of  sight).  If  N  =  1  then  the  missile 
is  burning  at  the  same  rate  as  the  LOS,  or  simply  homing  on  the  target.  If  N  <  1,  then  the 
missile  will  be  turning  slower  than  the  LOS,  thus  continually  failing  behind  the  target, 
making  an  intercept  impossible.  In  order  to  completely  understand  the.  physics  of 
proportional  navigation  guidance  it  is  necessary  to  analyze  pursuit  and  constant  bearing 
guidance. 

B.  PURSUIT  GUIDANCE 

For  pursuit  guidance,  the  missile  velocity  vector  is  always  directed  toward  the  target 
as  illustrated  by  Figure  2.1.  The  missile  is  then  constantly  heading  along  die  line  of  sight 
from  the  missile  to  the  target  and  its  path  describes  a  pursuit  path.  Given  that  the  rate  of  turn 
of  the  missile  is  always  equal  to  the  rate  of  turn  of  the  LOS,  “pure”  pursuit  (without  leading 
angle)  paths  are  highly  curved.  This  requires  the  missile  to  use  significant  acceleration. 


Since  the  signal  processing  is  limited  to  continuously  locating  the  target  and  changing  the 
missile  flight  path  angle,  the  on-board  avionics  are  relatively  simple.  As  will  be 
demonstrated  later,  this  kind  of  classical  guidance  law  is  a  special  case  of  PROPNAV  when 
the  effective  navigation  ratio  is  equal  to  1. 


Figure  2.2  shows  the  geometry  of  the  pursuit  guidance  law.  Vm  and  V t  are  respectively 
the  missile  and  target  velocities,  0m  and  8f  are  respectively  the  missile  and  target  flight 
path  angles,  ai  is  the  difference  between  the  LOS  angle  and  the  target  flight  path  angle,  and 


r  is  the  instantaneous  separation  between  missile  and  target.  Inertial  and  missile  translating 
coordinate  systems  are  also  shown  in  the  figure. 

The  velocity  of  the  target  with  respect  to  the  missile  is  given  by: 


v  =  Vt-Vm\  (Eq  2.2) 

v  =  re/.  +  r0c6.  (Eq2.3) 


Figure  2.2  Pursuit  Guidance  Geometry 


Writing  the  velocity  of  the  target  and  missile  in  terms  of  the  polar  unit  base  vectors  er  and 
ee,  we  get: 
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(Eq  2.4) 


Vt  =  Vlco%a[er-V!smatee, 

(Eq  2.4) 

=  Vm*r- 

(Eq  2.5) 

Equating  equations  2.2  and  2.3  and  uring  equations  2.4  and  2.5,  we  obtain: 

f  =  V,cos  arVm, 

(Eq  2.6) 

r®m  =  -v«sin“,- 

(Eq  2.7) 

From  Figure  2.2,  we  see  that: 

0m  =  a,  +  0,. 

(Eq  2.8) 

Considering  a  non  responsive  target: 

0m  =  a,. 

m  t 

(Eq  2,9) 

The  missile  acceleration  u  is  obtained  by  differentiating  equation  2.5: 


given  that  from  analytical  mechanics: 


(Eq  2.10) 


der  a  dQ 
dt  ~  e*di ' 


(Eq  2.11) 

Assuming  constant  speed  (magnitude  of  the  velocity  vector),  the  acceleration  command 
will  be  the  normal  component  of  the  acceleration  which  will  be  designated  um\ 

“m  -  Vm6m  =  VmX  =  Vmd,,  (Eq  2.(2) 

where  a,  is  a  time  function. 


C.  CONSTANT  BEARING  GUIDANCE 

The  accelerations  required  by  the  pursuit  guidance  law  can  be  reduced  by  aiming  the 
missile  ahead  of  the  target  by  using  a  lead  angle.  In  this  case  the  missile  traverses  a  straight 
line  to  a  collision  with  a  constant  speed  non  maneuvering  target,  as  shown  in  Figure  2.3. 
The  missile  converges  on  the  target  by  using  a  constant  LOS  angle  ( X  =  constant).  Since  the 
rate  of  change  of  the  LOS  angle  is  zero  throughout  the  flight,  the  lateral  accelerations  are 


zero.  If  the  target  maneuvers  evasively,  by  changing  its  velocity  vector  in  direction  and/or 
in  magnitude,  a  new  collision  course  must  be  computed  and  the  missile  flight  path  altered 
accordingly.  The  constant  bearing  geom^-ry  is  shown  in  Figure  2.4. 

We  wish  to  find  the  missile  control  input  necessary  to  responde  to  prescribed  target 
accelerations.  The  relative  velocity  is  given  by: 

V  =  V;- Vm  =  rer  +  rXee.  (Eq2.13) 

The  target  and  missile  absolute  velocities  are: 

V,  =  Vtcoscxter  -  V^sina^e;  (Eq2.14) 

Vm  =  Vmcosam^r-Vmsin(Xm^(i-  (Eq  2.15) 

Subtracting  equation  2.14  from  equation  2.15  and  equating  the  result  to  2.13,  we  find  that: 

r  =  V,  cos  a,-  V m  cos  am ;  (Eq2.16) 

rk  =  — Vj  sin  ot^  +  V  m  sin  ctm .  (Eq2.17) 

The  requirements  for  a  constant  bearing  guidance  are: 

X  =  0;  (Eq  2.18) 

r<  0.  (Eq  2.19) 

Using  equations  2.17  and  2.18,  we  get: 

since. 

sinocm  =  -y-Vn  (Eq  2.20) 

v  m 

From  equations  2.16  and  2.19: 

cos  a. 

cosam  >  ~y~~Vf  (Eq  2.21) 

in 
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Figure  2.3  Constant  Bearing  Guidance  Trajectory 


We  may  use  equation  2,20  to  obtain  a  expression  for  cosaf  by  squaring  the  equation  and 
using  the  fundamental  trigonometric  entity.  Doing  this,  we  get: 


cos  a 


t*. 


vi 


+  rcosu^;  - 


(Eq 


Figure  2.4  Constant  Bearing  Geometry 


Substituting  this  expression  into  equation  2.21,  we  obtain: 


cosoc„,  > 

m 


l 
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+  (cosam)2 

'  fJl ' 


(Eq  2.23) 


This  expression  is  satisfied  if: 


(Eq  2.24) 


V, 

-±-l<Q=>Vm>Vl 
* m 

Thus,  for  this  guidance  law  to  be  effective  the  missile  must  have  speed  advantage  relative 
to  the  target. 

The  missile  and  target  accelerations  can  he  computed  from  Figure  2.5,  which  expands 
the  acceieradon  in  terms  of  tangential  and  normal  components.  The  velocity  vector  of  a 
body  (missile  or  target)  is  described  by: 

V  =  Vx,  (Eq  2.25) 

where  x  represents  the  tangential  unit  vector  to  the  trajectories  represented  by  the  dashed 
lines  in  the  figure.  The  figure  shows  the  response  of  the  missile  to  a  evasive  target.  We  are 
interested  in  computing  the  relations  hip  between  the  missile  control  input  and  the  target 
maneuver. 
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The  acceleration  of  the  body  is  given  by: 


dV  ,1  dx 

u  =  ~  =  f(Vx)  =  Vx  +  V-j-. 
dt  dr  dt 


(Eq  2.26) 


For  small  values  of  Avy  it  approaches  the  magnitude,  of  Ax  and  the  direction  of  Ax  becomes 

dx 

perpendicular  to  the  direction  of  x .  It  follows  that  the  derivative  -  is  of  magnitude  1  and 

d\y 


perpendicular  to  x.  Then  this  derivative  is  the  unit  normal  vector  h .  The  time  derivative 


dx 

dt 


is  found  by  using  the  chain  rule,  as  follows: 

dx  dx  dy  . 

—  = - =  ny  =  nQ. 

dt  dyd't 


(Eq  2.27) 


Assuming  constant  speed,  the  missile  and  target  accelerations  are  always  in  the  direction  of 
the  respective  unit  normal  components  and  can  be  written  as: 

A m  =  (Eq  2.28) 

=  ntVfi(  (Eq  2.29) 


since, 


and  given  that  k  =  0, 


0m  =  a„  +  ?i, 
m  m 


0m  =  am, 

m  tn} 


u„  =  Vma  . 


m  F  m 


m 


(Eq  2.30) 

(Eq  2.31) 
(Eq  2.32) 


The  variable,  um  is  the  missile  acceleration  magnitude.  Differentiating  equation  2.20,  we 
find  am  and  then  um  in  terms  of  the  target  evasive  acceleration: 

(  cosa,(r)  v 

d»(,)  ■  <E«233) 
where  the  time  factor  is  included.  Additionally, 


LI 


a,(f) 

a. t(t)  -  0,(0  =  -y— 
v  t 


(Eq  2.34) 


where  u(  (t)  is  the  target  acceleration  magnitude.  Hence: 


(Eq  2.35) 


From  this  last  equation  and  equation  2.20,  we  conclude  that  the  LOS  will  maintain  its 


direction  in  space,  keeping  the  missile  on  a  collision  course  with  the  target  provided  that 
the  missile’s  and  target’s  kinematics  normal  to  the  LOS  behave  likewise.  Additionally, 
from  equation  2.21,  the  closing  velocity  (component  of  the  relative  velocity  along  the  LOS) 
must  be  positive.  Constant  bearing  guidance  requires  the  knowledge  of  the  heading  and 
velocity  of  the  target,  the  line  of  sight,  and  the  velocity  of  the  missile,  which  dictates  a  more 
complex  signal  processing  system  than  for  pursuit  guidance. 


D.  PROPORTIONAL  NAVIGATION  GUIDANCE 


1.  In  Search  Of  The  Proportional  Navigation  Concept 

Pursuit  guidance  tries  to  continuously  point  the  missile  to  the  target,  resulting  a 
highly  curved  path  and  very  large  accelerations.  The  guidance  law  is  only  interested  in  the 
present  position  of  the  target;  lacking  information  about  the  target  kinematics.  This  lack  of 
information  precludes  the  missile  from  building  a  lead  angle,  resulting  in  a  somewhat 
ineffective  guidance  law.  Constant  bearing  guidance  points  the  missile  to  the  future 
position  of  the  target,  resulting  in  a  straight  line  collision  path  with  a  non  maneuvering 
target.  Before  pointing  the  missile,  the  guidance  system  needs  to  know  the  heading  and 
velocity  of  the  target  to  compute  the  target’s  future  position.  So,  this  method  is  not 
practical,  especially  when  dealing  with  targets  with  evasive  capabilities. 

The  advantage  of  proportional  navigation  is  that  it  provides  a  practical  method  of 
approximating  a  constant  bearing  course  to  a  maneuvering  target.  PROPNAV  tries  to 
emulate  the  constant  bearing  guidance  command  by  using  LOS  rate  information  from  an 
on  -  board  electromagnetic  or  electro  -  optic  device. 
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A  missile  using  constant  bearing  guidance  only  needs  a  control  input  -when  it  is 
necessary  to  change  its  heatims?  at  the  begir  ning  of  the  flight  and  afterwards  if  the  target 
maneuvers.  The  form  of  this  command  signal  was  derived  and  is  repeated  here  for 
convenience: 


“m  =  Vma-mh- 


From  the  proportional  navigation  geometry  in  Figure  2.6: 


0m  =  a  +  X. 

m  m 


(Eq  2.36) 

(Eq  2.37) 


Taking  its  derivative: 

0m  =  am  +  iL  (Eq  2.38) 

The  acceleration  of  the  missile  (assuming  constant  speed)  is: 

fi*  -  W  -  ^  (Eq 2.39) 

Our  goal  is  to  emulate  equation  2.36  by  using  a  linear  transformation  between  the  LOS  rate 
X  and  the  missile’s  angle  rate  am.  Set,  for  example: 


i  1  • 

pj  —  l  nr 


Equation  2.39  becomes: 


-  V\*(a„  + 


1 


=  v. 


N  .  . 

nt  yY _  ni 


(Eq  2.40) 


(Eq  2.41) 


m  m  /v  —  1  m 

By  letting  N  be  large  this  equation  approaches  equation  2.36  for  a  constant  bearing  path 
(collision  course).  We  are  interested  in  the  relationship  between  the  LOS  rate  and  the  flight 
path  angle  rate.  Using  equations  2  38  and  2.40,  we  find  that: 


a  —  A.  .  1  —  /  A/  _  1  \  T  J  T  __  t  . 
M  t  a,  —  i/V  ii/vv/v  —  A  a, 

m  m  '  ' 


(Eq  2.42) 


um  =  NVmkn.  (Eq  2.43) 

Therefore,  PROPNAV  is  a  practical  guidance  law  that  emulates  constant  bearing  guidance 
by  issuing  control  commands  that  are  proportional  to  the  LOS  rate. 
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Pursuit  guidance  is  a  particular  case  of  proportional  navigation  when  N  =1 
(compare  equations  2.12  and  2.43).  As  we  have  seen,  constant  bearing  guidance  is  obtained 
by  letting  N  be  large  (theoretically  infinity)-  However,  large  gains  in  the  amplifiers  also 
cause  large  amplifications  of  noise;  therefore  N  is  usually  restricted  to  less  than  6. 
proportional  navigation  paths  are  less  curved  than  pursuit  paths,  but  more  curved  than 
constant  bearing  collisions.  PROPNAV  anticipates  the  future  position  of  the  target  without 
actually  computing  it.  Due  to  this  property  this  guidance  law  presents  a  higher  degree  of 
responsiveness  than  other  guidance  laws. 

2.  Proportional  Navigation  And  Zero  Effort  Miss 

In  Figure  2.6  the  missile,  with  velocity  magnitude!^,  is  heading  at  an  angle  of 

L  + HE  with  respect  to  the  line  of  sight.  The  angle  L  is  known  as  the  missile  lead  angle  and 
is  the  theoretically  correct  angle  for  the  missile  to  be  on  a  collision  triangle,  with  the  target. 
If  the  missile  is  launched  in  a  collision  triangle  with  a  non  evasive  target,  no  further 
accelerations  commands  will  be  required  to  hit  the  target.  The  angle  HE  is  known  as  the 
heading  error,  and  represents  the  initial  deviation  of  the  missile  from  the  collision  triangle. 

In  practice,  the  missile  is  usually  not  launched  exactly  in  a  collision  triangle,  since 
the  expected  intercept  point  is  not  known  precisely.  The  location  of  the  intercept  point  can 
only  be  approximated,  because  we  do  not  know  in  advance  what  the  target  will  do  in  the 
future.  In  fact,  that  is  why  a  guidance  system  is  required.  The  point  of  closest  approach  of 
the  missile  and  target  is  known  as  the  miss  distance.  Guidance  system  lags  or  subsystem 
dynamics  will  cause  miss  distance.  The  simplest  proportional  navigation  homing  loop  is 
shown  in  Figure  2.7  where  we  have  linearized  the  missile/target  engagement  by  using  the 
small  angles  approximation  (i.e.  we  assume  that  the  flight-path  angles  and  the  line  of  sight 
angle  are  small  in  order  to  linearize  the  engagement  geometry.  Then,  the  cosine  functions 
are  approximated  by  1  and  the  sine  and  tangent  functions  by  their  arguments).  In  a 
linearized  analysis,  the  range  equation  is  approximated  by  the  following  time  varying 
relationship: 
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(Eq  2.44) 


r  =  Vc(tF-[)  =  Vctgo 
where  Vc  is  the  closing  velocity,  tF  is  the  total  flight  time,  and  t  (time  to  go)  is  the  time 
until  the  end  of  the  flight.  As  shown  in  [Ref.  1]  the  miss  distance  will  always  be  zero  in  a 
zero-lag  proportional  navigation  homing  loop.  The  PROPNAV  guidance  law  used  in  the 
homing  loop  of  Figure  2.7,  and  also  the  most  used  in  the  literature,  is  not  the  one  derived 
in  equation  2.43,  but  the  following  one: 

um  =  NVcXnx  (Eq  2.45) 

where  fix  is  the  unit  vector  normal  to  the  LOS.  Then,  the  control  input  is  issued 

perpendicular  to  the  instantaneous  LOS.  It  can  be  easily  demonstrated  that  this  last 
expression  maintains  the  proportionality  between  the  missile  flight  path  angle  and  the 
angular  LOS  rate.  In  [Ref.  2],  Guelman  contrasted  “pure”  PROPNAV  (described  by 
equation  2.43,  wherein  command  accelerations  are  normal  to  the  missile  velocity  vector) 
and  “true”  PROPNAV  (described  by  equation  2.45,  wherein  command  accelerations  are 
normal  to  the  line  of  sight).  He  concluded  that  the  later  law  would  result  in  intercept  only 
if  the  initial  conditions  were  within  a  well-defined  subset  of  the  parameter  space.  In  the 
homing  loop  of  Figure  2.7,  the  seeker  provides  the  LOS  rate  by  taking  the  derivative  of  the 
geometric  LOS  angle.  The  noise  filter  processes  the  noisy  LOS  rate  measurements  to 
provide  an  estimate  of  the  LOS  rate.  The  guidance  command  is  generated  using  the  “true” 
PROPNAV  guidance  law.  The  guidance  system  must  cause  the  missile  to  maneuver,  by 
using  moving  control  surfaces.  The  seeker  and  the  guidance  system  dynamics  are  described 
by  differential  equations. 

The  presence  of  delays  in  the  homing  loop  creates  miss  distance.  In  the  presence 
of  guidance  system  dynamics,  the  heading  error  {HE)  and  target  maneuver  (target  evasive 
acceleration)  are  the  two  sources  of  miss  distance.  The  PROPNAV  guidance  law  can  be 
expressed  in  terms  of  the  zero  effort  miss.  The  zero  effort  miss  is  not  only  useful  in 
explaining  PROPNAV  but  is  also  useful  in  deriving  more  advanced  missile  control  laws. 
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Figure  2.6  Proportional  Navigation  Two  Dimensional  Engagement 

The  zero  effort  miss  is  the  distance  the  missile  would  miss  the  target  if  the  target 
continued  along  its  present  course  and  the  missile  made  no  further  corrective  maneuvers. 
Using  Figure  2.6: 

ZEMX  =  rx  +  vxtg0,  (Eq  2.46) 

ZEMy  =  ry  +  vytgo,  (Eq  2.47) 

where  ZEM  represents  the  zero  effort  miss,  r  is  the  missile/target  relative  distance,  and  v  is 
the  missile/target  relative  speed.  The  subscript  (x  or  y)  represent  the  projection  of  the 
respective  quantity  over  that  coordinate  axis. 


Id 


Figure  2.7  Zero  -  lag  Proportional  Navigation  Homing  Loop 
(Linearized  Engagement) 


The  ZEM  perpendicular  to  the  LOS,  is  given  by: 

ZEMplos  -  -ZEM  xsin\  +  ZEMy cosA. 
Using  equations  2.46  through  2.48,  we  obtain: 


t  ( r  v  —  r  v  ) 

zemplos  =  i° . yjL 


The  LOS  is: 


,  f  fM 

A  -  atan  ---  , 


(Eq  2.48) 


(Eq  2.49) 


(Eq  2.50) 


taking  its  derivative,  we  obtain: 


r2 


(Eq  2.51) 


* **■:  '-iV'Wv  '-xvv,. -'tT1 


Comparing  equations  2.49  and  2.51,  the  LOS  rate  may  be  expressed  in  terms  of  the 
component  of  the  zero  effort  miss  normal  to  the  LOS: 


ZEM, 


ZEMr 


V  t 

r  c  go 


(Eq  2.52) 


where  r  ■=  Vct  .  Then  the  PROPNAV  guidance  command  magnitude  can  It  expressed  in 


terms  of  the  ZEM  perpendicular  to  the  LOS: 


N  ZEM , 


um  — 

m 


(Eq  2.53) 


Thus,  we  conclude  that  the  PROPNAV  acceleration  command  that  is  perpendicular  to  the 
LOS  is  not  only  proportional  to  the  LOS  rate  and  closing  velocity  but  is  also  proportional 
to  the  zero  effort  miss  and  inversely  proportional  to  the  square  of  time  to  go.  The  efficiency 
of  PROPNAV  guidance  is  a  direct  consequence  of  this  dynamic  property.  This  is,  of  course, 
a  very  powerful  concept. 


IH.  DYNAMIC  IMAGE  PROCESSING 


A.  GENERAL 

A  missile  that  uses  a  TV  camera  and  a  seeker  (radar  or  IR),  or  instead,  two  TV  cameras 
is  considered.  A  setup  with  only  one  TV  camera  is  also  studied.  The  seeker  and  the  camera, 
or  the  two  cameras,  can  be  located  on  the  missile’s  nose  separated  by  a  transversal  distance 
d.  The  seeker  plus  the  single  camera  setup,  permits  the  missile  to  emulate  the  stereo  vision 
of  the  two  cameras  setup.  It  has  the  additional  advantage  of  tracking  the  target  at  the  early 
stages  of  the  engagement  using  solely  the  seeker’s  LOS  angle  information.  This  system 
permits  us  to  compute  the  3-D  target  motion  by  using  a  two  perspective  views  motion 
algorithm  and  the  target’s  spatial  direction  and  range  provided  by  the  seeker.  The  two 
cameras  setup  permits  us  to  use  image  plane  locations  in  two  views,  corresponding  to  a 

single  object  point  at  times  t,  and  t2,  to  determine  the  3-D  object  (target)  locations  X0  (t:) 

and  X0  (t2) .  The  one  camera  setup  also  permits  us  to  determine  the  motion  of  the  target,  as 
a  function  of  time.  This  is  done  by  using  a  two  perspective  views  motion  algorithm  and 
guessing  the  target’s  physical  dimensions  to  estimate  its  absolute  depth.  In  this  way,  we 
emulate  binocular  vision.  The  estimated  3-D  motion  of  the  target  and  the  image  sampling 
time  permit  us  to  estimate  the  target  velocity  and  acceleration  components  in  a  preselected 
3-D  rectangular  coordinate  system.  The  acceleration  information  can  subsequendy  be 
injected  into  the  control  algorithms,  which  will  be  developed  in  the  next  chapter,  to  improve 
the  miss  distance. 

1.  Scene  (3-D)  -  Image  (2-D)  Geometric  Considerations 

Mathematically,  we  can  express  the  transformation  of  object  point  locations  (3- 
D)  to  image  plane  locations  (2-D)  by  the  following  generally  noninvertible  geometric 
transformation: 

Xt(t)  =  g(X0  ('),...)•  (Eq  3.1) 
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The  modeling  of  the  imaging  process,  described  by  the  above  equation,  relates  object  points 

X0  (f)  in  the  3-D  scene  to  image  points  Xj  ( t)  in  the  image  plane.  The  function  g  depends 
on  the  imaging  geometry,  lens  model,  and  coordinate  system  choices. 

The  imaging  model  is  derived  by  considering  the  pinhole  camera  model  shown  in 

Figure  3.  1.  The  point  Xc  lies  over  the  camera‘s  optical  axis  at  a  distance /from  the  image 
plane.  The  figure  shows  two  distinct  coordinate  systems,  an  image  plane  coordinate  system 
and  a  global  coordinate  system.  Our  first  goal  is:  assuming  the  simplified  camera  model 
shown  in  Figure  3.1,  derive  the  transformation  described  by  equation  3.1  where  the  object 

point  X0  (r)  can  be  measured  from  either  coordinate  system. 
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The  object/image  relationship  defined  in  equation  3.1  is  defined  by  a  transformation 
matrix.  Independent  of  the  camera  model,  this  transformation  matrix  is  the  product  of  two 
matrices.  The  first  matrix  describes  the  object  -  image  coordinates  transformation,  and  is 
derived  by  assuming  that  the  3-D  object  coordinates  are  measured  relatively  to  the  image 
plane.  However,  if  the  object  coordinates  are  measured  relatively  to  the  global  coordinate 
system,  a  second  transformation  matrix  relating  the  two  coordinate  systems  have  to  be 
defined.  This  matrix  is  the  composite  of  the  relative  rotation  and  translation  between  the 
coordinate  systems.  It  describes  the  coordinates  transformation  between  the  two  coordinate 
systems. 

To  identify  the  transformation  defined  by  equation  3.1,  the  two  matrices  are 
derived  for  the  simplified  camera  model  of  Figure  3.1.  Monocular  vision  (only  one  camera) 
is  incapable  of  determining  absolute  depth;.  However,  any  imaged  point  is  constrained  to 

correspond  to  an  object  point  located  anywhere  on  the  3-D  line  segment  containing  Aj  ( t) ,  • 
Xb  and  X0(t) . 

Assuming  that  the  coordinate  systems  for  both  object  and  image  points,  are 
coincident  and  centered  in  the  image  plane,  the  above  colinear  points  are  related  by  : 

*(*<(')  -*c)  =  (Xc-X0(t)).  (Eq 3.2) 

Expa  iding  this  equation  yields: 

fc{[0y,.z^-[f0  0]7’}  =  [f0Q\r-\xoyoz^,  (Eq  3.3) 

where  the  superscript  T  is  the  transpose  operator.  The  time  index  t  has  been  dropped  to 
simplify  the  notation.  Equation  3.3  yields: 


1 

II 

'o  ^ 

II 

(Eq  3.4) 

y0  fy0 

y‘  "  k  ~  (/-*,)* 

(Eq  3.5) 
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The  minus  sign  in  the  second  expression  of  the  two  last  equations,  stands  for  the  image 
inversion  originated  by  the  back  -  projection  model  of  Figure  3.1.  The  matrix  representation 
of  the  nonlinear  equations  3.5  and  3.6  is: 


"o  /oo1 

Xj  =  MX0  =  0  o/oit  (Eq  3.7) 

-10  0/J 

The  last  equation  uses  homogeneous  coordinates  (a  technique  also  used  to  develop 
computer  graphics),  for  image  and  object  points.  The  homogeneous  coordinates  are  defined 
by  multiplying  the  physical  coordinates  by  an  arbitrary  constant  c  and  including  the 
constant  as  an  additional  element  of  the  vector; 

~  [cy,-  czl  c]  .  (Eq  3.8) 

and 


*o=\x0y0z0  l]r  (Eq  3.9) 

Note  that  in  equation  3.9  the  arbitrary  constant  is  equal  to  1.  The  object  -  image  point 
transformation  is  defined  by  equation  3.7,  where  it  is  implicitly  assumed  that  ,v(-  =  0 . 
Rewriting  equations  3.5  and  3.6,  we  conclude: 


(Eq  3.10) 


Fixing  the  image  plane  coordinates  y{  and  zt  the  above  equation  describes  the  3-D  line  over 
which  the  3  -D  object  is  located.  Therefore,  while  the  transformation  in  equation  3.7  is  not 
invertible,  choice  of  a  specific  image  point  constrains  corresponding  object  points  to  lie 
along  a  3-D  ray  (shown  in  Figure  3.1). 

If  the  object  points  are  measured  relatively  to  the  global  coordinate  system,  the 
matrix  relating  the  two  coordinate  systems  has  to  be  computed.  This  transformation  matrix 


is  the  product  of  a  succession  of  matrices.  Individually,  each  of  these  matrices  defines  a 
rotation  or  translation  of  the  image  plane  coordinate  system  relative  to  the  global  coordinate 
system.  The  succession  of  transformations  may  be  of  the  form; 

£0  =  (  T2  R2  R\  71  )X8o>  (Eq3.ll) 

where  X0  and  define  the  homogeneous  object  coordinates  in  the  image  plane  and  global 
coordinate  systems,  respectively.  Here  the  fust  transformation  is  the  translation  Tl 
followed  by  the  rotation  Rl,  etc.  The  composite  of  the  above  transformation  may  be  defined 

by: 

Hg->i  -  T2R2R1TI'-  (Eq3.12) 

Then,  the  general  relationship  between  object  points  measured  relatively  to  any  user 
selected  coordinate  system  and  the  image  plane  points  is: 

Xi  =  MHg^tXt.  (Eq  3,13) 

For  the  simple  case  of  only  a  translation  as  shown  in  Figure  3.1,  we  see  that  in  object 
coordinates: 

X0  =  Xl-[dxdyd]r.  (Eq  3.14) 

Homogeneous  coordinates  enable  us  to  represent  the  last  relationship  using  a  translation 
matrix: 


X0  =  H^iXl  = 


10  0-d, 
0  10  -dy 
0  0  1  -dz 
00  0  1 


(Eq  3.15) 


2.  Stereo  Vision  (2  Cameras) 

Monocular  vision  disables  depth  perception.  In  fact,  due  to  the  impossibility  of 
inverting  the  3  x  4  matrix  Q  =  obtained  in  the  last  section,  we  are  constrained  to 

the  determination  of  image  points  from  object  points.  However,  we  are  interested  in 
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determining  the  3-D  locations  (measured  relatively  to  a  global  coordinate  system).  One 
approach  to  soive  this  problem  is  to  use  more  than  one  camera.  One  of  our  proposals,  is  to 
use  two  cameras  in  the  missile’s  nose  separated  by  a  distance  d  emulating,  in  some  way, 
the  human  visual  system. 

Initially,  we  assume  the  simplified  two  dimensional  diagram  of  the  stereo  vision 
in  Figure  3.2.  The  scene  consists  of  a  2-D  surface.  As  shown  in  the  figure,  a  point  on  this 
surface  is  projected  onto  the  two  image  planes  (IP1  and  IP2).  In  general  the  two  centers  of 
projection  differ  in  length  (/j  and  /2).  It  is  assumed  that  the  user  selected  global  coordinate 
system,  to  measure  the  object  coordinates,  is  coincident  and  centered  in  the  image  plane 
IP1.  The  coordinate  xi  is  shown  in  the  figure,  the  coordinate  yi  is  perpendicular  to  xi  and 
in  the  plane  of  the  page.  The  object  point  may  be  determined  using  the  two  projected  points, 
one  in  each  camera. 

The  relationship  between  the  homogeneous  coordinates  of  the  object  point, 
measured  relatively  to  the  image  plane  IP1,  and  the  homogeneous  coordinates  of  the 
corresponding  image  point  is: 


p- 

~i  o  o' 

*0 

cxil 
_  c . 

— 

^0 

1 

(Eq  3.16) 

The  object  coordinates  measured  from  IP2  are  related  to  the  object  coordinates  measured 
from  IP1  by  the  following  relationship: 


or,  in  homogeneous  coordinates: 
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(Eq  3.17) 


(Eq  3.18) 
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Figure  3.2  Two  Dimensional  Stereo  Vision 


Then  image  points  in  IP2,  denoted  xi2,  may  be  related  to  object  points  measured  relatively 
to  the  global  coordinate  system  by: 


■■  1 

1  0  0 

'l  0  d 

*0 

cxi2 

= 

0  1 

0  1  0 
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L  h  . 
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(Eq  3.19) 


Using  equations  3.16  and  3.19  the  following  relationships  in  object  coordinates  are 


obtained: 
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xn(fi  -y0) 
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(Eq  3.20) 
(Eq  3.21) 


Equating  these  two  equation  the  object’s  depth  yQ  may  be  found: 
>0  ”  fi*n  ~fixi 


(Eq  3.22) 


The  object  point  x0  may  be  found,  from  either  equation  3.20  or  equation  3.21.  Hence,  using 

two  image  planes  permits  us  to  determine  the  object  point  depth  yu  from  its  corresponding 

image  points.  This  was  proved  for  a  2-D  surface.  Next  we  are  going  to  see  how  to  do  it  in 
a  3-D  environment. 

Equation  3.13,  defmes  the  relationship  between  the  scene  three  dimensional 
points  and  the  correspondent  two  dimensional  image  points  by  using  a  matrix  of 
transformation: 


Q  =  (Eq  3.23) 

Hg  ;  is  the  coordinate  systems  transformation  matrix  which  depends  on  the  rotations  and/ 
or  translations  of  the  image  plane  coordinate  system  relative  to  the  user  selected  global 
coordinate  system.  M  is  the  object  -  image  transformation  matrix,  which  is  a  function  of 
the  imaging  geometry  and  lens  model.  The  Q  matrix  is  a  non  -  invertible  3x4  matrix 
(assuming  homogeneous  coordinates)  and  may  be  genetically  represented  by: 


/-» 

14  = 


^11  112113  1 14 
121  122  123  12A 


111  1 32  1 33  1 34j 


(Eq  3.24) 


The  missile  must  have  sufficient  processing  capability  to  find  this  matrix  in  real  time.  The 
object  -  image  transformation  matrix  M  is  generally  invariant  (however  zooming  the  scene, 
for  example,  changes  its  value).  The  coordinate  systems  transformation  matrix  Hg  f  must 
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be  dynamically  updated  as  the  missile/target  engagement  proceeds.  The  selected  global 
coordinate  system  for  missile  guidance  simulation  purposes  is  a  ground  coordinate  system 
which  will  be  presented  in  the  next  chapter.  For  full  dimension  (3  -D)  stereo  vision,  the  two 
cameras  arrangement  may  be  described  using  homogeneous  coordinates  as: 

Xic  =  QCX0<  (Eq  3.25) 

where  the  index  c  =  1,2  refers  to  the  cameras.  Since  each  of  these  two  matrix  equations 
(one  for  each  sensor)  represents  two  equations  in  physical  coordinates,  we  obtain  four 
equations  and  three  unknowns  by  using  the  two  cameras  stereo  arrangement.  The  matrix 
equation  3.25  can  be  explicitiy  written  for  each  sensor  as: 


and. 
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(Eq  3.26) 
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(Eq  3.27) 


The  index  of  each  matrix  element  is  composed  by  three  numbers,  the  first  is  the  camera 
number  and  the  next  two  represent  the  element  position  into  the  matrix.  Each  of  these  two 
matrix  equations  generates  two  equations  in  physical  coordinates.  To  find  these  equations, 
the  arbitrary  constants  (cL  and  c2)  have  to  be  calculated.  Then,  each  of  the  constants  is 
substituted  into  the  two  remaining  matrix  equations.  Finally,  regrouping  terms  as 
coefficients  of  xg0,  yg0  and  zg0  ,  a  set  of  four  equations  is  obtained.  Performing  this  procedure 
to  obtain  the  first  physical  equation  from  the  matrix  equation  3.26,  we  get: 

cl  =  4131*0  +  4132^  +  4133zo  +  4i34-  (Eq  3.28) 
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substituting  this  expression  into  and  regrouping  terms,  we  obtain: 


(Eq  3.29) 

((<7ni“4t3i.v/iK+  (</i2i~^i32.va)>’o+  ’ 

The  set  of  four  equations  and  three  unknowns  is  written  compactly  in  matrix  notation  as: 

PXi  =  F  (Eq  3.30) 

or 


4lll  “4l31.vil  4121  ~4i32>,il  ^113  “9i33>’il 

4 134  ^/l  ~  4 U4 

^121  “</l31zil  4l22'“4l32z/l  4i23“4l33zii 

xi  = 

4 134  zil  ~4l24 

4211  “  4231  y(2  4212  ~  4232^/2  4213  ~  4233>12 

4234  y-il  ~  4214 

4221  “4231z/2  4222  “  4232z/2  4223  “  4233z»2_ 

4234  2 i2  ~  4224 

Equation  3.31  may  be  solved  using  least  square  techniques  by  forming  the  pscudoinverse 
of  P  denoted  /*+ .  Hence: 


Pf0  =  F  =t>  PTPtl  =  PrF  =>  f0  -  (PTP)  ~lPTF  =>xl  =  P\F.  (Eq  3.32) 

This  equation  yields  the  mean  square  estimate  for  the  object  point  a0.  Alternatively.  Xxu 
may  be  found  by  using  three  of  the  four  equations,  assuming  that  the  three  equations  are 
linearly  independent. 

In  this  exposition,  we  have  assumed  that  the  necessary  image  plane  point 
correspondences  have  been  determined.  It  is  important  to  say  that  this  is  the  most  difficult 
problem  in  the  development  of  a  stereo  vision  algorithm.  Techniques  to  solve  this  problem 
are  presented  in  [Ref.  3].  [Ref.4]  presents  an  algorithm  to  match  stereo  images. 
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B.  ESTIMATING  3-D  MOTION  PARAMETERS  OF  A  RIGID  BODY  FROM 
TWO  CONSECUTIVE  IMAGE  FRAMES 


I.  General 

In  section  A  of  this  chapter,  we  have  shown  that  using  two  cameras,  we  are  able 
to  find  the  object  (target)  3-D  position  relative  to  a  user  selected  global  coordinate  system. 
As  intermediate  steps,  it  is  necessary  to  establish  feature  correspondences  between  selected 
points  in  the  two  stereo  images  (static  stereo).  It  is  also  necessary  to  establish  feature 
correspondences  for  all  pairs  of  consecutive  image  frames  in  each  camera’s  image 
sequence.  The  targets  that  we  are  interested  in  are  mainly  airplanes.  Therefore,  we  may  use 
as  points  for  feature  correspondences  the  tips  of  the  wings,  nose,  stabilizers  and  rudder. 

Estimation  of  the  3-D  target  acceleration  components  may  be  divided  in  three 
steps.  In  ihc  first  step,  the  target  estimated  points  at  ^  and  f2  are  used  to  estimate  its  3-D 


velocity  components.  In  the  second  step,  the  3-D  velocity  components  of  the  target  are 
computed  using  the  target’s  estimated  points  at  t2  and  /3.  Finally,  the  third  step  estimates 
the  target’s  3-D  acceleration  components  by  identifying  the  time  change  of  the  target’s 
velocity  for  each  of  the  three  velocity  components.  Alternatively,  we  may  use  Kalman 
filtering  theory  to  estimate  the  3-D  target  acceleration  components. 

A  different  approach  for  estimating  the  3-D,  time  dependent,  target  motion  is  now 
presented.  The  formal  definition  of  a  rigid  3-D  object  is  one  for  which  the  3-D  distances 
between  any  pair  of  points  on  the  object  do  not  change  with  time  that  is,  for  all  pair  of  points 
on  the  3-D  object: 


Y 


"  cmn< 


Vf,  Vm, 


(Eq  3.33) 


where  cmn  are  constants.  The  assumption  of  an  rigid,  or  nondeformable  target  is  reasonable 
and  creates  additional  constraints  for  motion  estimation.  Rigidity  constrains  the  motion  of 
individual  object  points  to  be  strongly  coupled,  although  the  need  for  point  correspondence 
maintains.  Then,  the  3-D  translation  of  the  target  may  be  determined  by  estimating  the 
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translation  parameters  of  a  single  point  object.  The  basis  of  estimating  the  target’s  motion 
using  this  approach,  is  that  the  3-D  motion  of  a  rigid  body  can  be  described  by  a  3-D 
translation  vector  and  three  rotation  angles  chosen  with  respect  to  a  user  selected  coordinate 
system,  Then,  six  parameters  completely  define  the  target’s  motion.  Formulating  the 
rotations  using  three  rotation  matrices  (Rq,  Ra  and  R  ),  the  target’s  motion  is  described 

by: 


X0(t2)  =  RX0(tl)  +f,  (Eq  3.34) 

where  R  is  the  overall  rotation  matrix: 


R  -  /?0/?a/?p, 


(Eq  3.35) 


and  T  is  the  translation  matrix.  Equation  3,34  may  be  represented  in  homogeneous 
coordinates  as: 
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0  0  0 
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x0(h). 

lj 


(Eq  3.36) 


2.  Monocular  Motion  Estimation  Using  Two  Perspective  Views 

Our  goal  is  to  compute  dynamically  the  rotation  (JR)  and  translation  (T)  matrices 
from  point  (or  feature)  correspondences  between  two  perspective  view’s.  We  can  divide  the 
process  of  estimating  the  three  -  dimensional  motion  of  the  target  from  image  sequences  in 
three  steps.  The  first  step  is  to  establish  feature  correspondences  between  two  consecutive 
image  frames.  Correspondences  between  features  may  be  established  through  matching  or 
inter  >  frame  hacking.  [Ref. 4]  develops  a  two  -  view/stereo  matcher  that  computes 
displacement  fields  from  two  images.  The  second  step  is  to  estimate  the  motion  parameters 
(R  and  T  matrices).  The  third  step  is  to  estimate  the  3-D  motion  of  the  target  using  equation 
3.34.  There  are  a  number  of  papers  in  the  digital  image  processing  and  computer  vision 
literature  dealing  with  estimation  of  the  motion  parameters  of  the  3-D  motion  of  a  rigid 
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body  from  two  consecutive  image  frames.  Weng,  Huang  and  Ahaja  [Ref.  5],  propose  an 
algorithm  that  given  3  point  correspondences  solves  for  a  intermediate  matrix  called  the 
essential  parameter  matrix  (E).  Then  the  Rotation  matrix  (/?)  and  the  translational  direction 

t 

(the  unit  vector  — - )  are  obtained  from  E.  The  magnitude  of  the  translational  vector  (II  T'i ) 

* !  i  il 

and  the  absolute  depths  of  the  object  points  ( xok  and  x‘ ok  where  xok  is  the  absolute  depth 
of  the  object’s  k  feature  point  and  x' ok  is  the  absolute  depth  of  the  object’s  k  feature  point 
after  being  rotated  and  translated)  cannot  be  determined  by  monocular  vision.  The 

xok  x'ok 

algorithm  also  solves  for  the  relative  depths  and  ).  The  algorithm  is  unable  to 

estimate  the  target’s  position  because  it  is  not  possible  to  calculate  the  absolute  depth.  This 
agrees  with  intuition,  due  to  the  lack  of  invertibility  of  the  3-D  to  2-D  image  transformation. 
To  overcome  this  problem,  we  propose  to  estimate  the  absolute  depth  by  correlating  the. 
dimensions  of  the  target  over  the  image  plane  with  the  guessed  physical  dimensions  of  the 
target.  A  relatively  easy  trigonometric  approach  permits  us  to  estimate  the  target’s  depth 
given  the  target’s  physical  dimeni.  ;ons.  Another  approach  for  emulating  stereo  vision  is  to 
use  range  and  directional  spatial  information  from  the  seeker.  This  information  is  combined 
with  the  monocular  vision  equations  in  order  to  estimate  the  target’s  3-D  motion. 

In  conclusion,  monocular  vision  may  be  applied  to  estimate  the  motion  of  the 
target  if  additional  information  about  the  target  is  available  (or  guessed).  The  target’s 
acceleration  components  may  now  be  estimated  and  injected  into  the  missile’s  control 
algorithm. 

3.  Stereo  Motion  Estimation  Using  Two  Perspective  Views 

As  we  have  seen,  two  or  more  spatially  distributed  sensors  enables  determination 
of  the  3-D  target  motion.  Here,  our  goal  is  to  determine,  not  only  the  target’s  motion  but 
also  the  rotation  and  translation  matrices  that  describe  the  motion.  Given  two  consecutive 


time  samples  or  “frames”  in  each  sensor,  as  well  as  image  point  correspondences,  we  may 
write: 


x!(0  =  Gj-;*o(0, 

X?(t)  =  Q^jXoiO, 


(Eq  3.37) 
(Eq  3.38) 


x!('  +  ts) 


=  Q'.Xoit  +  rj  =  q 
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\X0(t) ,  (Eq  3.39) 


0  0  0  |1 


fi(t  +  Ts)  =  Q2g^iX0(t  +  Ts)  =  Q2g_l 


R 


0  0  0|1 


x0(t) .  (Eq  3.40) 


The  matrix  t  is  the,  already  known,  3-D  to  2-D  transformation  matrix  (the  superscript 
refers  to  the  sensor).  Ts  is  the  time  between  two  consecutive  images.  The  system  is 
represented  in  homogeneous  coordinates.  Assuming  n  point  for  point  correspondences,  a 
total  of  8/j  equations  in  physical  coordinates  is  obtained.  However,  the  number  of 


unknowns  is  12  +  3n  (9  elements  of  R,  3  elements  of  f  and  3  elements  for  each  X0).  This 
yields  the  constraint  on  the  number  of  point  for  point  correspondences: 

12  +  3  n£Sn.  (Eq3.41) 

Since  n  must  be  an  integer,  n  £  3 .  Thus,  three  corresponding  image  points  from  two  views 
in  two  frames  are  sufficient  to  determine  both  the  motion  parameters  and  the  3-D  location 
of  the  object  points. 

The  three  target  acceleration  components  may  be  computed  by  taking  the  second 
derivative  following  the  filtering  of  the  target’s  motion  data.  Alternatively,  the  target’s 
motion  may  be  processed  by  Kalman  filters  to  estimate  its  acceleration  components. 


This  additional  information  about  the  target’s  behavior  may  be  used  to  improve 
the  missile  guidance  towards  the  target.  In  order  to  effectively  use  this  information,  we  have 
first  of  all  to  determine  control  laws  that  can  use  and  produce  better  results  if  this 
information  is  available.  This  will  be  stressed  in  the  next  chapter. 
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rv.  SIMULATION  DEVELOPEMENT 


A.  CONTROL  ALGORITHMS  DEVELOPMENT 

Chapter  n  introduced  proportional  navigation  guidance.  In  this  section  we  derive  more 
advance  guidance  laws.  Contrary  to  PRQPNAV,  these  guidance  laws  use  the  estimated 
acceleration  of  the  target  as  a  additional  input  to  the  homing  loop.  As  will  be  seen  in 
Chapter  V,  the  advanced  guidance  laws  relay,  the  interceptor  acceleration  requirements  and, 
in  general,  yield  smaller  miss  distances. 


I.  Augmented  Proportional  Navigation 

Proportional  navigation  issues  control  commands  that  are  proportional  to  the 
predicted  zero  effort  miss  normal  to  die  line  of  sight  ( ZEMPL0S ).  That  is,  the  missile 
guidance  system  tries  to  minimize  the  final  miss  distance  between  the  target  and  the  missile 
by  issuing  acceleration  commands  that  arc  proportional  to  the  miss  distance,  that  would 
result  if  the  missile  made  no  further  corrective  acceleration  and  the  target  did  not  maneuver. 
Therefore,  if  the  target  maneuvers  evasively  it  generates  additional  miss  distance  that  is  not 
accounted  for  in  the  PROPNAV  guidance  law.  Augmented  proportional  navigation  also 
issues  guidance  commands  that  are  proportional  to  the  predicted  miss  distance.  However, 
for  augmented  PROPNAV,  the  miss  distance  is  estimated  by  taking  into  account  the 
maneuver  of  the  target  (target  acceleration).  The  augmented  PROPNAV  target's 
acceleration  dependent  term  will  be  calculated.  This  term  is  injected  into  the  homing  loop 
to  enhance  the  guidance  performance.  In  the  following  analysis,  we  follow  the 
nomenclature  of  Chapter  II  and  the  geometry  of  Figure  2.6. 

The  x  component  of  the  miss  distance,  for  an  evasive  target,  is  computed  as 
follows  (the  y  component  is  computed  similarly): 

^  (MO)  =  MO.  (Eq 4.1) 

where  t'  represents  time.  Then, 
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(Eq  4.2) 


h 

ZEMx(l)  -  rx[lm['  =  ^(0  +jvx(t‘)d[‘. 

t 

Where,  Z EMx(t)  is  the  x  component  of  the  zero  effort  miss,  predicted  at  time  t'  -  t  and 
rx  ( t )  is  the  present  missile  -  target  relative  distance  along  the  x  axis.  But, 

<(0  =  ~  (MO),  (Eq 4.3) 

where  a*  (/')  is  the  x  component  of  the  target  acceleration.  Then, 

V«(f)  r 

|  dvx(t')  =  ja*(r)dt",  (Eq 4.4) 

V,  (0  t 

where  r"  is  the  variable  of  integration.  Hence: 

v 

vx{t')  =  vx(t)  +  ja*(t")dt".  (Eq  4.5) 

t 

Substituting  this  equation  into  equation  4.2,  we  get  the  expression  for  the  predicted  x 
component  of  the  ZEM  at  time  f. 

hr 

ZEMx{t)  =  rx\(m(f  =  rx(t)  +  vx(t) /g0  +  JJ< (f) dt" dt" .  (Eq4.6) 

1 1 

Where,  tgo  =  tF  - t  is  the  time  to  go.  The  y  component  of  the  ZEM  (t)  is  obtained  using 
the  same  reasoning,  and  is: 

hr 

ZEMy(t)  =  ry\tmtf  =  ry(t)  +  vy(t)  tg0  + jjay  (t")  dt" dt"  .  (Eq4.7) 

t  i 

r  r 

The  two  interior  integrals  ( t” )  dt "  and  j ayt  (t” )  dt "  are  time  functions;  call: 

f  r 

r 

kxV)  =  jax,(t")dt",  (Eq  4.8) 

t 

and 
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ky(n  =  ft  in  dr. 


Then,  equations  4.6  and  4.7  may  be  expressed  as: 


(Eq  4.9) 


(Eq  4.10) 


ZEMx(t)  =  rx\[mtf  -  M‘')  +vx(t)tg0  +  ^kx(tH)dr 


=  rx(t)  +Vx(t)  tgo  +  hx(tgo) 


(Eq  4. 1 1) 


ZEMy(t)  =  ry\(mt'  =  ry(t)  +  vy(t)tg0  + ft(t")dtu  -  ry(t)  +vy(t)tg0  +  hy(tg0)\ 


where  hx{tg0)  and  hy  (tg0)  are  tg0  and  maneuver  dependent  functions.  The  component  of  the 
ZEM  that  is  perpendicular  to  the  LOS,  ZEMPL0S,  is: 

ZEMPL0S(t)  mZEMy(t)Ga*(\(t))-ZEMx{t)&in(\(t)).  (Eq4.12) 

rJt) 

Substituting  4.10  and  4.11  into  this  equation  and  setting  cos(A(r))  =  and 


sin(X(0)  = 


ry(t) 

TU) 


,  we  obtain: 


ZEMptosiO  - 


(Eq  4.13) 


r  (r)  L  (0 

(r  (r)  +v  (r)t  +  h  (t  ))— ^ — (r  (/) +v  (t)  t  +h  (t  ))-->-r- 

\  y  \  i  yv/go  yVgo,,j-/A  vxv/  xv/go  X  80  f(/) 


From  equation  2.52: 


MO 


rx(t)vy(t)  -  ry(t)vx(t ) 
r2  (t) 


(Eq4.14) 


Then  equation  4.13,  which  takes  into  account  the  target’s  acceleration  to  estimate  the  miss 
distance,  may  be  reduced  to: 
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rM)  r  (t) 

ZEMpL0S(t)  =  tgor(t)XU)  +  ( hy(tg0 )  --(7y  -hx(tgo)  ~-)  •  (Hq4.15) 

The  PROPNAV  guidance  command  may  be  expressed  in  terms  of  the  ZEM  as  (see  derivation 
in  Chapter  0): 


um(0 


NVc(t)  ZEM,  >s(/) 
”  ~rU)  lg0 


(Eq  4.16) 


where  ZEM  PLOS(t)  was  then  derived  for  a  nonresponsive  target,  If  the  target  maneuvers  the 
zero  effort  miss  is  augmented  by  an  additional  term,  on  the  right  hand  side  of  equation  4.15. 
Therefore,  a.  perfectly  plausible  guidance  law,  in  the  presence  of  target  maneuver,  would  be: 

(Eq  4.17) 
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m'APN 


(0  =  NVc{t)h(t)  + 


NVAt) 

vrtTt~  (hy  ^ W  cos  ( X  ( /) )  -  hx  ( t  )  sin  ( X  ( t) ) )  • 

'  lgo 


This  guidance  law  is  PROPNAV  with  an  extra  term  that  accounts  for  the  maneuver  of  the 
target.  The  equation  was  derived  for  a  nonlinearized  geometry.  The  impossibility  of  knowing, 
a  priori,  the  future  target  maneuver,  precludes  the  calculation  of  h y(tgo) ,  hx(t  )  and  tg0. 
However,  if  the  target  desires  to  inflict  the  most  miss  distance  it  must  maneuver  at  a  small  time 
to  go.  Also,  considering  the  time  constants  associated  with  the  target’s  maneuver,  we  propose 
to  approximate  the  time  dependent  target  acceleration  at(t')  by  a  constant  target 
acceleration^,  (t) )  i.e.assume  at  (/’)  =  a,  (  t) .  Then,  the  control  law  may  be  stated  as: 

(Eq  4. 18) 


u 


m  APN 


(!) 


=  NVc(t)X(  t) 


r\t)t 


s° 


or. 


•  Wet™ 

=  NV(X  +  — (aJcos(X)  -a'sin(X)) ;  (Eq4.19) 


'  m  APN 
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where  the  time  factor  was  dropped.  The  extra  term,  present  in  augmented  PROPNAV 
(expression  in  parenthesis  on  the  right  hand  side  of  equation  4.19),  is  proportional  to  the 
component  of  the  target  acceleration  normal  to  the  LOS. 

Linearization  of  the  nonlinear  missile  -  target  geometry  is  shown  in  [Ref.  1]  to  be 
an  accurate  approximation  to  the  actual  geometry.  Then,  assuming  a  linearized  geometry, 
equation  4.19  may  be  further  reduced  to: 

um  =  NV^  +  ^ay,  (Eq  4.20) 

where  we  have  considered  that  the  LOS  angle  is  small.  A  zero  -  lag  augmented  proportional 
navigation  homing  loop  assuming  linearized  geometry  is  shown  in  Figure  4. 1 . 


Figure  4.1  Zero  -  lag  Augmented  Proportional  Navigation  Homing  Loop 


38 


The  additional  target  maneuver  term  required  by  die  guidance  law,  appears  as  a 
feedforward  term  in  the  homing  loop  block  diagram. 

2.  Optimal  Intercept  Guidance 

The  missile  -  target  engagement  scenario  may  be  described  in  state  space 
representation  by  the  following  linear  system: 

HO  -f[x(t),u( r),r]  =  Ax(t)  +  Bu(t);  (Eq 4.21) 

x  is  the  n  -  dimensional  state  vector  describing  the  relative  movement  between  the  missile 
and  the  target  and  also,  the  dynamics  of  the  guidance  system.  The  variable  u  is  the  m  - 
dimensional  missile’s  control  input  vector.  We  seek  to  find  a  guidance  law  that  is  a  function 
of  the  system  states.  There  is  an  infinite  number  of  possible  guidance  laws.  Thus,  it  is 
necessary  to  state  in  mathematical  terms  what  the  guidance  law  should  do.  Certainly  we 
wish  to  design  a  terminal  controller  that  would  bring  certain  components  of  ,v  ( tF)  to  zero, 
using  “acceptable”  levels  of  control.  One  way  to  do  this  is  to  minimize  a  performance  index 
made  up  of  a  quadratic  form  in  the  control: 

r,  tF 

J  =  JZ,[jc(0,w  (/),*]<//  =  2\u2(t)dt,  (Eq  4.22} 

0  0 

subject  to  the  terminal  constraint: 

xt(tF)  =0,  i  =  0,1 . p,  (Eq  4.23) 

and  the  constraints: 

no  =/[*(;),  u(r),r]  =Ax(0+Bu(0,  (Eq4.24) 

x  (0) ,  given.  (Eq4.25) 

In  equation  4.23,  p  ^  n . 

The  miss  distance  will  always  be  zero  in  a  zero  -  lag  PROPNAV  navigation 
homing  loop.  Guidance  system  lags  or  subsystem  dynamics  will  cause  miss  distance. 
Optimal  guidance  eliminates  miss  distance  by  canceling  out  the  guidance  system  dynamics. 
In  this  way  the  optimal  guidance  law  attempts  to  make  the  real  world  guidance  system 
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appear  to  be  a  “perfect”  (zero  -  lag)  guidance  system.  To  find  the  optimal  control  vector, 
u  (/) ,  that  brings  the  system  from  a  initial  state  .v  (0)  to  a  terminal  state  x(tF)  (where 
some  of  its  components  are  zero),  we  can  use  the  method  of  the  Lagrange  multipliers.  Then 
the  constraints  (4.23)  and  (4.24)  may  be  adjoined  to  the  performance  function  (4.22)  by 

using  the  multipliers  s  -  (e, . e  ,  0p  +  1 ....  Gn) T  and  X  =  (X( . Xfl) 7  as  follows: 

i  '' 

J  =  eTx(tF)  +if{n2(r)  +  Xr(r)  [Ax(t)  +Bu(t)  -x(t)  ] }  cit.  (Eq  4.26) 

0 

The  Hamiltonian  is  defmed  as  follows: 

H[x(t),u(t)tt]  =  L[x(t),u(t),t] +\Tf[x(t),u(t),t] .  (Eq 4.27) 

Integrating  the  last  term  on  the  light  hand  side  of  equation  4.26  yields: 

(Eq  4.28) 
h 

J  =  zTx(tF)  -  kr(tF)x(tF)  +  Xr(0)x(0)  +J  {H[x(t),u{t)tt]  +X(t)x{r)  }  dt- 


Considering  the  variation  in  J  due  to  variations  in  the  control  vector  u  (0 ,  we  get: 

(Eq  4.29) 


ty 


dJ  =  [(er-A,7,)4te]l.^+  [Xrdx]/-0+j|  (^  +  kT{t))dx  +  °£du\dt 

0 

In  order  to  make  the  variations  in  J  due  to  variations  in  u(t)  independent  from  the 
variations  in  x(t)  produced  by  die  variations  in  u  (r) ,  we  choose  the  influence  functions 
X  (t)  to  cause  the  coefficients  of  dx  to  vanish: 

j . 


dH 


V(i) 

ox  ox  ox 


Then: 


X(r)  =  -ATX(t), 


with  boundary  conditions: 


X(tr-)  =  e. 


(Eq  4.30) 

(Eq  4.31) 

(Eq  4.32) 
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Using  these  results,  equation  4.29  becomes: 

h 

dJ  =  Xr(0)dx(0)  +j| ~dudt.  (Eq  4.33) 

o 

Hence,  XT  (0)  is  the  gradient  of  J  with  respect  to  variations  in  the  initial  conditions,  while 
holding  u  ( t)  constant  and  satisfying  the  constraints  of  the  problem.  For  an  extremum,  dJ 
must  be  zero  for  arbitrary  du  ( t ) .  This  can  only  happen  if: 

l  =  0,  0  (Eq  4.34) 


or: 

ur  +  XrB  =  0.  (Eq  4.35) 

Then,  we  may  determine  the  control  vector  u(t) ,  as: 

u(0  =  -BrX(S) .  (Eq  4.36) 

Substituting  equation  4.36  into  equation  4.24  and  repeating  equation  4.31,  the  following 
two  -  point  boundary  value  problem  is  obtained: 


fii 

1a 

Ja 

!_o 

-bbt \  [«" 

-Ar\[K 

■ 

(Eq  4.37) 

The  2 n  boundary  conditions  are: 

x  (0) , 

given, 

(Eq4.38) 

MW 

=  0, 

i  =  1,  . 

•  ■.p: 

(Eq  4.39) 

w  - 

0, 

i  -  p+  1 . n. 

(Eq  4.40) 

The  n  boundary  conditions  4.39  and  4.40  may  be  replaced  by  the  boundary  condition  4,32, 
which  may  be  rewritten  as: 

Mrf)  "  ei:  i  =  (Eq  4.41) 

\(tF)  =  0;  <=p  + (Eq  4.42) 
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The  two  -  point  boundary  value  problem  4.37,  4.38,  4.41  and  4.42  may  be  solved  by  the 
sweep  method  [Ref.  6].  The  sweep  method  seeks  to  find  solutions  of  the  form: 

k(t)  =  W(t)xU)  +Y(t)z;  (Eq  4.43) 

z  =:  U(t)x(t)  +V(t)E.  (Eq 4.44) 

where  W (t) ,  Y (t) ,  U  (t)  and  V (t)  are  time  dependent  matrices, 

-  [.v. . .vJ  .  and  e;  =  [e, . ej  =  [k . k J  ,  .  Therefore,  we  want 

1  r  l  m  1  / ’  1  l'  t  m  lf 

to  find  solutions  for  the  influence  functions  k  ( t)  that  are  function  of  the  state  vector  .v  ( /  ) 
and  the  final  value  of  the  influence  functions,  or  equivalently,  of  the  specified  final  states 
z.  Since  equations  4.43  and  4.44  must  be  valid  at  t  =  tF\ 

W(tF)  =  0,  (Eq  4.45) 

=0,  (Eq  4.46) 


=  - 
i  0 


( n-p )  xp 


■■  nxp 


(Eq  4.47) 


U(tF) 


‘ P*P 


0 


p  X  (n-p) 


■  i 


! 

‘  /;  x  n 


(Eq  4.48) 


where  I  is  the  identity  matrix  and  0  is  a  zero  matrix  with  the  specified  dimensions. 
Substituting  equation  4.43  into  4.37  and  treating  e  as  a  constant  vector,  we  get: 

Wx  +  Wx  +  re  =  -4 T  ( VT.c  +  fe) .  (Eq  4.49) 

Substituting  x  from  equation  4.37  into  the  last  equation,  and  again  using  equation  4.43  to 
eliminate  A,  we  obtain: 

(W  +  WA+ATW-WBBTW)x+  (ATY  +  Y-WBBTY)e  =  0.  (Eq  4.50) 


This  expression  must  be  true  for  any  .v  and  e,  so  the  coefficients  of  x  and  £  must  vanish: 
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(Eq  4.51) 


W+  WA  +ATW-WBBTW  -  0;  W(tF)=  0. 

and 

Y  +  ArY-WBBrY  =  0:  Y(ty)=-  -  (Eq4.52) 

®ti  -/>  Xp 

n  x  /) 

Next,  we  differentiate  equation  4.44  with  respect  to  time,  treating  a  and  :  as  constant 
vectors: 

Ux  +  Ux  +  Ve  =  0.  (Eq  4.53) 

Substituting  x  from  equation  4.37,  and  using  equation  4,43  to  eliminate  A.,  we  obtain: 

( 0  +  UA  -  (JBBrW)  x  +  (V  -  UBBrY)  e  =  0.  (Eq4.54) 

This  last  expression  must  be  true  foi  any  .*  and  e,  so  the  coefficients  of  ,v  and  e  must 
vanish: 

U  +•  UA  -  UBBTW  -  0,  (Eq  4.55) 

V-UBBrY  =  0.  (Eq  4.56) 

From  equations  4.52  and  4.55  and  the  boundary  conditions  4.47  and  4.48,  we  conclude  that: 

£/(/)=  Yr(t).  (Eq  4.57) 

Then  equation  4.56  may  be  rewritten  as: 

V  =  YTBBTY\  V(tF)  =  0.  (Eq  4.58) 

The  Ricatti  equations  4.51,  4,52  and  4,58  may  be  integrated  backwards  from  the  final 
conditions  to  yield  iV(r),  Y  (r)  and  V(t) .  The  equation  4.44  is  solved  for  c  to  yield: 

£  =  \V(t)r[[z-YT(t)x(t)].  (Eq  4.59) 
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Our  goal  is,  to  find  the  influence  functions  X(t)  in  order  to  find  the  optimal 
control  vector  using  equation  4.36.  Equation  4.59  may  be  substituted  ifuo  equation  4.43  to 
find  X  (/) ; 


X(t)  -  (W-nrlK7').v(r)  +YV'z. 
However,  the  Ricatti  equation  4.51  has  as  solution: 

W{t)  =  0. 

Hence,  equations  4.52  and  4.58  become  simply: 


Y  +  ArY  =  0; 


and 


r 

I 


Ipxp 

^ll-pxp 


(Eq  4.60) 

(Eq  4.61) 

(Eq  4.62) 


h' 

V{t)  =  -j(YTBBTY)dl.  (Eq  4.63) 

i 

Combining  equations  4.36, 4.60  and  4.61  yields; 

u  (;)  =  -BrX  (t)  =  (~BrYV~l )  [z  -  Yrx(t)  ] .  (Eq  4.64) 

The  final  condition  that  we  are  interested  on  is  z  -  [xp  ....  xp]  f'm  f  =  [0 . 0]/. 

Hence  the  expression  4.64  may  be  reduced  to: 

u(0  =  BrYV~lYrx(i),  (Eq  4.65) 

where  Y  and  V  are  computed  from  equations  4.62  and  4,63,  respectively. 


Now  that  we  have  derived  the  terminal  controller  optimal  feedback  guidance  law, 
we  proceed  to  derive  the  continuous  feedback  law  described  by  equation  4.65  for  a  single 
lag  guidance  system.  The  single  -  lag  guidance  system  is  mathematically  described  in  the 


Laplace  domain  as: 
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1 

ur? 


(Eq  4.66) 


where  um  is  the  missile’s  acceleration,  um  is  the  command  acceleration,  and  T  is  the 

effective  guidance  system  time  constant. 

The  relative  motion  between  the  target  and  missile  is  considered  with  the 
linearized  (small  angles  approximation)  intercept  geometry  shown  in  Figure  4.2.  The 
assumption  of  small  angles  (flight  path  angles  0m,  0f  and  LOS  angle  \)  permits  us  to 

express  the  equations  of  motion  in  terms  of  state  variables  normal  to  the  reference  intercept 
course.  The  single  -  lag  guidance  model  shown  in  Figure  4.3  integrates  the  missile  -  target 
relative  motion  of  Figure  4.2  with  the  dynamics  of  the  guidance  system.  The  diagram  of 
blocks  of  Figure  4.3  can  be  expressed  in  state  space  form  as: 
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(Eq  4.67) 
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equations  4.62  and  4.63,  respectively.  The  solution  for  Y  is  obtained  from: 


Y  +  A  Y  =  0,  Y(tF)  = 


l‘i 

!o , 

j°! 

!oi 


(Eq  4.69) 


;o  o  o  o  i 

|i  o  o  o 
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Y  -  0,  Y  (tF)  = 


!  1 

|° 

10 

i 

:0 


Y  is  obtained  by  integrating  equation  4.70  backwards,  and  is: 


(Eq  4.70) 


r  = 


lF  1 


go 

to 


| -T2 

[  Jj:  t 

e  T  +  —  1 

L  7  Jj 

L 

T 

(Eq  4.71) 


L 

The  matrix  V  (in  this  case,  V  is  a  scalar  since  we  are  specifying  only  one  terminal 
constraint,  namely  zero  miss  distance:  y  =  0),  is  computed  using  equation  4.63: 

(Eq  4.72) 

h  h  /  V  -  2 


V (/)  »  -jYTBBTYdt  =  (I  ;  J  ^  r  +  -  1^ 


dt\  where  tg0-tF-t 


After  some  cumbersome  computations,  we  find: 


V  (t)  =  V  ("  ^  +  3*2  -  3 k)  +  12 ke~k  +  3e~2k  -  3;  where  k  = 


(Eq  4.73) 
T 


Then  using  equation  4.65  we  obtain  the  optimal  feedback  control  law: 

u(t)  =  ~  [y  +yt  +G.5d'i2  +  (~amT2)  [e~krk-  1)] ,  (Eq4.74) 


go 


where  N  = 


6k2  (e~k  +  k  -  1) 


zxii  k=  r  ■  •  (Eq4-75) 


-k  n  -2  k ’ 


T 


2k 3  -  6k2  +  6k  +  3  -  12 ke~K  -  3e 

It  is  desirable  to  express  the  state  variables  y  and  y  in  terms  of  the  line  of  sight  rate  A. 
Assuming  small  angles,  we  find  from  Figure  4.2  that: 
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(Eq  4.76) 


(*  =  {  = 
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X  - 
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c  50 


r  Vc(tF-t) 

hence,  tiic  optimal  control  law  can  be  expressed  as: 

u(t)  =  NVck+^{e-k  +  k-l)am  +  "ar 

h. 


(Eq  4.77) 


This  equation  is  a  biased  proportional  navigation  guidance  law  where  a  time  varying 
navigation  gain  N  and  an  acceleration  feedback  path  provide  compensation  for  the  missile 
time  lag.  The  acceleration  command  is  issued  normal  to  the  LOS. 


B.  TRIDIMENSIONAL  MISSILE/TARGET  ENGAGEMENT 

In  this  section  we  are  going  to  model  the  missiie/target  tridimensional  engagement 
scenario.  Three  guidance  laws,  namely:  PROFNAV,  augmented  PROPNAV  and  optimal 
guidance  will  be  used  and  tested  in  missile  guidance.  The  engagement  for  the  two  later 
guidance  laws  will  be  modeled,  by  assuming  the  presence  of  a  seeker  and  a  camera  aboard 
the  missile  to  extract  the  target’s  3-D  acceleration.  The  three  dimensional  MATLAB 
programs  are  presented  in  Appendices  A  through  C.  The  simulation  results  are  presented 
in  the  next  chapter. 


1.  3-D  Missile  /Target  Geometry 

The  tridimensional  scenario  is  developed  in  spherical  coordinates  by  defining  two 
perpendicular  planes  in  pitch  and  yaw,  as  illustrated  by  Figure  4.4. 

In  Figure  4.4  r  is  the  reladve  distance  between  missile  and  target;  A,  and  A,V(iW 

are  the  line  of  sight  angles  over  the  pitch  and  yaw  planes,  respectively,  and  may  be 
expressed  as: 


\itch  =  atan 


(■7—7  ) 


(Vr-.V/J2. 


(Eq  4.78) 
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KaW  --  atanS  TT^T) 

L  xm> 


(Eq  4.79) 


The  coordinate  system  shown  in  Figure  4.4  translates  with  the  missile.  To  track 
the  missile  target  tridimensional  positions  we  define  a  ground  based  coordinate  system 
shown  in  Figure  4.5. 


In  Figure  4.5: 
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m_  pitch 


=  atan 


V  >  v2  +  y2 


jy 

m' 


(Eq  4.80) 
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pitch 


=  atari 


■rf+yt 


(Eq  4.81) 


Figure  4.5  Ground  Coordinate  System 

®<i4-82> 

\.,m  =  atan(',|)  (Eq  4  83) 

The  missile  is  controlled  in  3-D  space  by  issuing  guidance  commands  in  two 
orthogonal  planes  (pitch  and  yaw),  upi!ch  and  uvaw.  The  magnitude  of  these  commands 
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depends  on  the  selected  guidance  law  and  their  tridimensional  direction  is  perpendicular  to 
the  LOS  defined  in  the  two  planes  (see  Figure  4.4).  The  LOS  in  pitch  is  defined  by  the 
imaginary  3-D  line  from  the  missile  to  the  target  in  the  pitch  plane.  The  LOS  in  yaw  is 
defined  by  the  imaginary  2-D  line  from  the  missile  to  the  projection  of  the  target  over  the 
yaw  plane.  The  yaw  plane  is  simply  the  horizontal  xy  plane.  The  pitch  plane  is  the  vertical 
plane  normal  to  the  horizontal  plane  and  rotated  by  the  yaw  angle  kvaw. 

The  guidance  laws  under  study  can  be  expressed  as  a  function  of  the  classical 
proportional  navigation  guidance  law  plus  a  term  that  may  depend,  among  other  variables, 
on  the  target  and  missile  accelerations.  Hence,  each  guidance  law  (classical  proportional 
navigation,  augmented  proportional  navigation  and  optimal  guidance)  can,  in  general,  be 
expressed,  as: 

um(t)  =  NVcX+f(amtaf,tg0,T),  (Eq4.84) 

where  um  (r)  is  the  guidance  command  that  is  issued  perpendicular  to  either  the  LOS  in 
yaw  or  the  LOS  in  pitch.  N  is  constant  for  PROPNAV  and  augmented  PROPNAV.  For 
optimal  guidance,  N  is  function  of  the  time  to  go  and  the  effective  time  constant  of  the 
guidance  system.  The  closing  speed  Vc,  is  the  relative  speed  between  the  target  and  the 

missile  along  either  the  LOS  in  yaw  or  the  LOS  in  pitch.  The  LOS  rate  k  may  be  the  LOS 
rate  in  either  the  pitch  or  the  yaw  planes.  The  term  f{atn>  ar  t  ,  T)  may  be  function  of  both 
the  missile’s  acceleration  am  and  the  target’s  acceleration  ar  the  time  to  go  i^0  and  the 
guidance  system’s  effective  time  constant  T.  In  order  to  generate  the  pilch  and  yaw 
guidance  commands  'upitcfl  (r)  and  u  (t) ,  it  is  fust  necessary  to  explain  how  to  obtain 
the  variables  that  they  depend  on. 

2.  Seeker  Head  Modeling 

The  seeker  is  able  to  detect,  acquire  and  track  by  sensing  and  processing  the 
radiation  or  reflection  of  energy  by  the  target.  The  seeker  is  normally  located  in  the 


51 


missile’s  nose  and  mounted  on  a  gimballed  platform  which  maintains  the  target  within  the 
field  of  view  by  rotating  the  platform. 

The  control  torque  to  the  seeker  may  be  described  by  the  following  equation. 

T  =  /p.  (Eq  4.85) 

where  T  is  the  applied  torque,  I  is  the  seeker’s  moment  of  inertia  and  j3  is  the  seeker’s 
angular-  acceleration.  The  seeker’s  dynamics  is  modeled  by  the  following  second  order 
differential  equation: 

13  =  ]  =  -c-j  (P-A.)  ~c2P.  (Eq  4,86) 

where  the  coefficients  cq  and  c2  are  determined  by  the  seeker's  time  constant  (t  .)  and 

damping  ratio.  Taking  the  Laplace  transform  of  equation  4.86,  assuming  zero  initial 
conditions,  we  obtain  the  filter’s  transfer  function  that  represents  the  relationship  between 
the  LOS  angle  input  A.  (.v)  and  the  seeker  head  angle,  output  |3  (s) : 


PCs) 

Ms) 


S  +  CoS  +  C 


(Eq  4.87) 


Assuming  a  damping  ratio  of  one,  the  transfer  function  in  equation  4.87  may  be  rewritten 
as: 
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(Eq  4.88) 


Choosing  xyi  =  0.1  sec  (which  is  a  good  approximation  of  a  real  world  system),  the 
constants  c,  and  c~,  may  be  obtained: 


c,  =  (-1)  =  100, 


\vi 


■2  =  2(-i-)  = 

\sk 


20. 


(Eq  4.89) 


(Eq  4.90) 
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Given  that  we  are  interested  in  the  3-D  missile/target  engagement,  the  seeker  must  provide 
line  of  sight  rate  information  in  both  planes.  Hence: 


$  pitch  __  20 

^ pitch  s^  +  20s  +100 

=  20 

^vaiv  +  20s  + 100 


(Eq  4.91) 
(Eq  4.92) 


where  &pitch  and  (3vaiv  are  the  seeker’s  pitch  and  yaw  angles,  respectively.  Figures  4.6  and 
4.7  depict  the  pitch  and  yaw  signal  flow  graphs  of  tire  seeker. 


From  these  diagrams,  the  continuous-time  state  equations  of  the  form: 


*sk  ^■skxsk'^'  ^ sk^sk’ 


(Eq  4.93) 


may  be  easily  obtained.  Selecting  the  state  vector  to  be: 


Figure  4.7  Seeker  Head  Flow  Graph  (Y aw) 


and  the  seeker  head  input  as: 


equation  4,93  becomes: 
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(Eq  4.95) 


(Eq  4.96) 


The  variables  $pilch  and  $yaw  are  estimates  of  the  LOS  angle  rates  kpitch  and 

Xyaw,  and  are  available  from  the  second  and  forth  states  of  xsk,  respectively.  The  estimates 
of  the  LOS  rate  in  pitch  and  yav/  permits  us  to  determine  the  missile  command  inputs  upiich 


3.  Guidance  System 

In  this  work,  the  guidance  system  dynamics  are  modeled  as  a  single  lag  as  seen  in 
equation  4.66,  This  equation  is  repeated  here. 


amJs) 

“«  (s) 


\  4-  Ts 


(Eq  4.97) 


For  the  3-D  missile/target  engagement  the  guidance  system  generates  missile  commands  in 
both  planes,  pitch  and  yaw.  Hence: 


am.  pitch  (■*')  _  1 

“pitchTsT  "  TTf  s ' 


am .  yaw 
^ yaw 


1  4  7V 


(Eq  4.98) 

(liq  4.99) 


where  am  pUch  and  am  yaw  are  the  pitch  and  yaw  missile’s  accelerations;  upitch  and  uyaw 

ai'e  the  pitch  and  yaw  missile’s  acceleration  commands.  Figures  4.8  and  4.9,  show  the  pitch 
and  yaw  signal  flow  graphs,  for  the  missile  guidance  system.  We  chose  the  guidance  system 
time  constant  T  to  be  1.0  second. 


5j 


Figure  4.8  Guidance  System  Signal  Flow  Graph  (Pitch) 


Figure  4.9  Guidance  System  Signal  Flow  Graph  (Yaw) 


From  these  diagrams,  the  state  equation  is  easily  obtained.  Defining  the  guidance  system 
state  vector  as: 

_  am_  pitch 
'  8*  ~  a 

h  m.  yaw 

and  the  guidance  system  input  as: 


(Eq  4.100) 
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(Eq  4.101) 


u  =  \UPitch\ 

V  , 


the  state  equation  becomes: 

_  -l  0! 


A‘.VV  = 


ltva\v 1 


,  +;iol 


0  -11  *5  '0  11 


(Eq  4,102) 


4.  Missile  And  Target  Kinematics 


The  missile  is  controlled  in  three  dimensional  space  by  generating  acceleration 
commands  in  two  orthogonal  planes.  These  planes  are  the  pitch  and  yaw  planes.  The 
acceleration  commands  in  pitch  and  yaw  are  issued  perpendicular  to  the  respective  lines  of 
sight.  The  magnitude  of  the  acceleration  commands  depends  on  the  selected  guidance  law. 
From  equation  4.84,  the  pitch  and  yaw  missile  acceleration  commands  may  be  expressed, 
in  its  general  form,  as: 


^ pitch  c _  pitch^pitch  pitch  Llt'  ^  4.103) 

liyuw  ~  yaw^'yuw^'fyuw^m' ■  (Eq  4.104) 

Vc  pilcit  and  Vc  yaw  are  the  relative  speeds,  between  the  target  and  the  missile  along  the 
pitch  and  yaw  line  of  sights.  The  functions  fpitch  («m,  a,,  tg0,  T)  and  fyav(am ,  cir  tg0,  '/') 

are  the  augmented  PROPNAV  or  optimal  guidance  extra  terms  in  pitch  and  yaw, 
respectively. 

In  order  to  track  the  missile’s  3-D  coordinates  (jcm,  ym,  zm) ,  the  missile  command 

accelerations  in  pitch  and  yaw,  are  broken  down  into  cartesian  coordinate  system 
components. 


Figure  4.10,  shows  the  decomposition  of  the  pitch  acceleration  command  in  its 
components.  From  this  figure  the  following  relationships  are  derived: 

amx.  puch  =  “  («m.  pi,chs'm\i,d)  C0SSnw’  ^  4' l05) 

amx.  pitch  =  -tem.pitchS^puch)  ^yav’  ^  4-106) 
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(Eq  4.107) 


a mi.  pitch  “  am.  pitch^os^pitch' 


Figure  4.11  shows  the  decomposition  of  the  yaw  plane  acceleration  command  in 


to  its  x  and  y  components.  From  the  figure  the  following  relationships  are  obtained: 

(Eq  4.108) 

(Eq  4.109) 


^mx.yaw  ~  ~~am_  vmvs^^ycjiv’ 


amy_yuw  (im.  yawc(i^^yaw' 
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To  find  the  overall  missile’s  acceleration  components  aiong  the  three  cartesian  axis,  we  use 
the  results  of  the  equations  4.105  through  4.109: 


'tnx  ~  ^mx.  pitch'*' ^mx.  yaw 

(Eq  4.110) 

*my  ~  umy.  pilch  Umy.  yaw' 

(Eq  4.111) 

umz  ~  amz.  pitch' 

(Eq  4. 112) 

59 

The  missile’s  tridimensional  movement  is  determined  by  these  three  acceleration 
components.  Defining  the  missile  state  vector  as: 


mi 


Xm  = 


•  mi 

:Vmi 

■zm\ 

::mj 


(Eq  4.1 13) 


and  the  input  as  the  missile’s  acceleration  command: 


n  1 
lwmxl 

am  ~  j amy 

L 


(Eq  4.114) 


the  missile  state  equation  is: 
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(Eq  4.115) 


As  we  have  seen  in  Chapter  HI,  the  missile  when  equipped  with  a  camera  and  a 
seeker  is  able  to  estimate  the  target’s  3-D  acceleration  components.  This  information  may 
be  used  to  improve  the  missile  guidance  towards  the  target.  Defining  the  target  state  vector 
as: 


Gc" 

n 


X,  = 


!3V| 

.  i  ■ 

i?'i 

rn 

! 


(Eq  4.116) 


and  the  target’s  acceleration  as: 
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(Eq  4.1 17) 


.  !X{ 

at  =  \aty\ - 

I  l 


the  target  state  equation  is: 


X,  =  i 


0  1  0  0  0  01  OOOl 
OOOOOO!  .  1  0  Oj 

ooo  i  o oi Y  +  ;oooi^. 
ooooooi' '  -0  1  0  •’ 
ooooo  i|  |o  o  o 

ooooooi  !o  0  1 


(Eq  4.118) 


the  tridimensional  target  acceleration  at  may  be  e  mated  using  the  missile’s  image 
processing  capabilities; 

5.  Pitch  And  Yaw  Closing  Velocities.  Determination  Of  Time  To  Go 

Figure  4. 12  shows  the  decomposition  of  the  missile  and  target  absolute  velocities. 
From  this  figure: 

I  T/  ^ 


V„  - 


m 


mx\ 


"Oj  • 


(Eq  4.1 19) 


V,  = 


V  l 

i  mz  i 


’  lx 

vv\ 

lv«i 


(Eq  4.120) 


Figures  4. 13  and  4. 14  show  the  projection  of  the  missile's  velocity  vector  over  the  pitch  and 
yaw  planes. 

These  figures  permit  us  to  compute  the  missile’s  and  target’s  velocity 
components,  over  the  pitch  and  yaw  planes: 

Vm.pitch  =  Vm™sUm_yaw~Kaw)<  (Eq  4.121) 


=  Vmcosr, 


(Eq  4.122) 
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(Eq  4.123) 


ma 


.pitch  ^rCOS  (”^f_  vavv  \vaw^ 
^.ymv  =  ^C0SY,.  ver' 


(Eq  4.124) 


Figure  4.12  Missile  And  Target  Velocity  Components 


The  pitch  closing  velocity  V.  pitch  is  found  by  projecting  the  missile’s  and  target’s  pitch 
plane  velocities  along  the  LOS  in  pitch  (see  Figure  4,4),  Then: 

(Eq  4.125) 

Vc.  pitch  ~  K m .  pitchC0S  (  ^ pitch  ~  Ym_  pilc^  ' ' ^t.  pitch C0S  (  ^ pitch  ~  \ 
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Similarly,  the  yaw  plane  closing  velocity,  is  obtained  as: 


(Eq  4.126) 


V'„ 


c_  yaw  ^ m_  (Y m.  yaw  ^ yaw ^  ytm>cos  ( Yt  yaw  ^yaw) 


Figure  4,13  Missile’s  Pitch  And  Yaw  Velocity  Components 


The  time  to  go  until  interception  may  be  computed  from: 

<*„=  \r—-  (Eq  4.127) 

v  c.  pitch 

where  r  is  the  missile/target  relative  distance  along  the  pitch  LOS  (see  Figure  4.4). 
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Figure  4.14  Target’s  Pitch  And  Yaw  Velocity  Components 

6.  Proportional  Navigation 

The  PROPNAV  3-D  simulation,  in  Appendix  A,  uses  missile  commands  in  pitch 
and  yaw  of  the  form: 


''pitch  pitch^pitch' 

(Eq  4.128) 

U  =  NV  k 

yaw  r  c.  yaw-yaw 

(Eq  4.129) 

A /  i  o  a  /'*rvr»cfarit 
WilotW  l  ’  iO  (A  VOUtiUUit! 

7.  Augmented  Proportional  Navigation 

The  augmented  PROPNAV  3-D  simulation,  in  Appendix  B,  uses  missile 
commands  in  pitch  and  yaw  of  the  form: 


“yaw  =  ^V'c_  yaw\aw  +  0.5 Nat_  yaw  (Eq  4.13 1) 

where  at  piuh  and  at  yaw  are  the  components  of  the  target  acceleration  normal  to  the 
pitch  and  yaw  line  of  sights  and  N  is  a  constant. 


8.  Optimal  Guidance 

The  optimal  guidance  3-D  simulation,  in  Appendix  C,  uses  missile  commands  in 
pitch  and  yaw  of  the  form: 


N 


u pitch  NVc.  pitch^" pitch  ,  2  ^  1)  am_  pitch  2  Ci<-  pitch’  (Eq  4.132; 

V-  =  NVC (Eq  4.133) 

K 

where  k  and  N  are  given  by  equation  4.75. 

9.  Discrete-Time  Simulation  Using  State  Space  Methods 
The  general  continuous-time  state  equations  are: 

i  U)  =  Ax(t)  +Du{t),  (Eq  4.134) 

v  (0  =  Cx(t)  +  Du  (r)  (Eq  4.135) 

where  x  ( t )  is  the  state  vector  and  y(t)  is  the  output  vector.  This  system  is  simulated  by 
iterating  the  discrete-time  state  equations: 

.r(«+l)  =  0x(n)  +  H/  (n) ,  (Eq4.136) 

>■(/!)  =  Cx(fi)  +Du(n).  (Eq  4.137) 

where: 


AT 

=  e  Ts  =  sampling  time: 

T, 

A  =  j/'fl  dt. 

0 


(Eq  4.138) 
(Eq  4.139) 
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The  missile/target  engagement  scenario  is  simulated  using  the  MATLAB 
software  package.  The  discrete-time  state  equations  used  in  the  simulation  are  defined  for 
the  seeker,  guidance  system,  missile,  and  target  dynamics: 


xsk  ( n  +  1 )  =  <&skxsk  ( n )  +  T %kusk  (n) ;  (Eq  4. 140) 

xgs(n  +  1)  =  <&g!iXgS{n)  +  r gsugs(n)  ;  (Eq  4,141) 

xm  (n  +  1)  =  <Um.vm  (n)  +  rmum  (n):  (Eq  4.142) 

x,  (n  +  1)  -  <!>,.*,  (n)  +  T tut(n) .  (Eq  4.143) 
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V.  SIMULATION  RESULTS 


A.  GENERAL 

This  chapter  presents  the  results  of  the  computer  simulations  for  a  missile  equipped 
with  a  radar  and  a  camera.  The  missile/target  3-D  engagement  was  simulated  using  three 
control  laws; 

1 .  Proportional  navigation, 

2.  Augmented  proportional  navigation, 

3.  Optimal  guidance. 

The  simulation  was  conducted  for  two  different  target  maneuvers: 

1 .  A  3-D  constant  target  acceleration, 

2.  A  3-D  varying  target  acceleration. 

The  following  assumptions  are  made  throughout: 

1.  We  assume  that  the  simulation’s  initial  conditions  are  defined  when  the  missile 
enters  into  the  terminal  phase  of  the  flight  (about  10  seconds  before  impact). 

2.  The  PROPNAV  and  APROPNAV  effective  navigation  constant  is  3  (A  =  3). 

3.  The  missile  is  limited  to  25  g’s  accelerations  in  pitch  and  yaw, 

4.  The  instantaneous  target  acceleration  is  available  from  previous  image 
processing.  No  delays  are  assumed  in  this  process. 

5.  The  missile  and  target  speeds  are  limited  to  3500  and  2000  feet/sec, 
respectively, 

6.  The  missile  is  in  a  collision  triangle  with  the  target  on  entering  into  the  terminal 
phase  of  flight. 

7.  The  target  may  start  its  evasion  maneuver  at  any  time  (initial  time), 

8.  The  target  is  limited  to  a  maximum  of  12  g’s, 

9.  The  acceleration  due  to  gravity  is  ignored, 

10.  The  sampling  time  is  0.01  seconds. 


B.  ENGAGEMENT  SCENARIOS.  RESULTS 


1.  Scenario  #1  (Constant  Target  Acceleration) 

The  initial  missile/target  geometry  (when  the  missile  enters  into  the  terminal 
phase  of  flight)  is  shown  in  Figure  5.1. 


fhe  engagement  initial  conditions  are  (distances  are  in  feet  and  speeds  nn  feet/ 


sec): 
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(Eq  5.2) 


The  target’s  evasive  maneuver  is  constant  (the  accelerations  are  in  feet/secA2) 


|3x« 

?t 

=  lx* 

V 

[?  x  8, 

(Eq  5.3) 


where  the  acceleration  of  gravity  is:  g  =  32.2  feet/sec2.  Figures  5.2  to  5.22  display  the 
results  of  the  three  dimensional  simulation  for  the  constant  target  evasive  maneuver. 
Figures  5.2  to  5.8  relate  to  the  proportional  navigation  (PROPNAV)  guidance  law.  Figures 
5.9  to  5. 15  relate  to  the  augmented  proportional  (APROPNAV)  guidance  law.  Figures  5.16 
to  5.22  relate  to  the  optimal  guidance  law.  Figures  5.4  to  5.8,  5.11  to  5.15  and  5.18  to  5.22 
display  the  results  assuming  that  the  target  starts  its  maneuver  6  seconds  after  the  missile 
entered  into  the  terminal  phase  of  flight. 
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Figure  5.4  Missile  Acceleration  Magnitude  (PROPNAV) 


Figure  5.5  Target  Acceleration  Magnitude  (PROPNAV) 
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Figure  5.7  Missile  Yaw  Acceleration  (PROPNAV) 
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Figure  5.9  Miss  Distance  vs.  Time  To  Go  (APROPNA  V) 


Figure  5.10  Miss  Distance  vs.  Initial  Time  (APROPNAV) 


MISSILE  ACCELERATION  MAGNITUDE  vs  TIME,  APROPNAV 


Figure  5.11  Missile  Acceleration  Magnitude  (APROPNAV) 


TARGET  ACCELERATION  MAGNITUDE  vs  TIME,  APROPNAV 
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Figure  5.12  Targe*  Acceleration  Magnitude  (APROPNAV) 
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Figure  5.13  Missile  Pitch  Acceleration  (APROPNAY) 


Figure  5.14  Missile  Yaw  Acceleration  (APROPNAY) 
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Figure  5.15  Missile  To  Target  Range  (APROPNAV) 


Figure  5.16  Miss  Distance  vs.  Time  To  Go  (OPTIMAL) 


Figure  5.17  Miss  Distance  vs.  Initial  Time  (OPTIMAL) 
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Figure  5.18  Missile  Acceleration  Magnitude  (OPTIMAL) 


Figure  5.19  Target  Acceleration  Magnitude  (OPTIMAL) 
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ure  S.20  Missile  Pitch  Acceleration  (OPTIMAL) 


Figure  5.21  Missile  Yaw  Acceleration  (OPTIMAL) 
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Figure  5.22  Missile  To  Target  Range  (OPTIMAL) 
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2.  Scenario  #2  (Varying  Target  Acceleration) 

The  initial  missile/target  geometry  and  the  initial  conditions  are  the  same  as  used 
for  tne  constant  acceleradon  scenario.  Namely,  the  target  geometry  is  shown  in  Figure  5.1 
and  the  initial  conditions  are  defined  by  equations  5.4  and  5.5.  The  target’s  evasive 
maneuver  may  start  at  any  time  to  go.  The  target  performs  the  following  3-D  maneuver  (the 
accelerations  are  in  feet/secA2): 

(Eq  5.4) 

:  3  x  c  x  sin  (  y,  )  j 

|  '■•t.yaw'  | 

I  4  x  8  x  cos  (Y,  vaJ 

|3  X  £  X  COS  (y,  pitch) 

where  yl  yaw  and  yt  iKh  are  the  target’s  yaw  and  pitch  flight  path  angles,  respectively. 

Figures  5.23  to  5.43  display  the  results  of  the  three  dimensional  simulation  for  this  type  of 
evasive  maneuver.  Figures  5.23  to  5.29  relate  to  the  proportional  navigation  (PROPNAV) 
guidance  law.  Figures  5.30  to  5.36  relate  to  the  augmented  proportional  (APROPNAV) 
guidance  law.  Figures  5.37  to  5.43  relate  to  die  optimal  guidance  law.  Figures  5.25  to  5.29, 
5.32  to  5.36  and  S.39  to  5,43  display  resuits,  assuming  that  the  target  starts  its  maneuver  6 
seconds  after  the  missile  entered  into  the  terminal  phase  of  flight. 
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Figure  5.23  Miss  Distance  vs.  Time  To  Go  (PROPNAV) 


Figure  5.24  Miss  Distance  vs.  Initiai  Time  (PROPNAV) 
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Figure  5.25  Missile  Acceleration  Magnitude  (PROPNAV) 


Figure  5.26  Target  Acceleration  Magnitude  (PROPNAV) 
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Figure  5.28  Missile  Yaw  Acceleration  (PROPNAV) 
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Figure  5.29  Missile  To  Target  Range  (PRQPNAV) 
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Figure  5.30  Miss  Distance  vs.  Time  To  Go  (APROPNAV) 


Figure  5.31  Miss  D'Stance  vs.  Initial  Time  (APROPNAV) 
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MISSILE  ACCELERATION  MAGNITUDE  vs  TIME.  APROPNAV 


Figure  5.32  Missile  Acceleration  Magnitude  (APROPNAV) 


Figure  5.33  Target  Acceleration  Magnitude  (APROPNAV) 
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Figure  5.34  Missile  Pitch  Acceleration  (APROPNAV) 


Figure  5.35  Missile  Yaw  Acceleration  (APROPNAV) 
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ure  5.39  Missile  Acceleration  Magnitude  (OPTIMAL) 
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Figure  S.41  Missile  Pitch  Acceleration  (OPTIMAL) 
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VI.  CONCLUSIONS  AND  RECOMEND ATIONS 


A. 


CONCLUSIONS 

1.  Scenario  #1  (Constant  Target  Acceleration) 


Figure  5.44  Miss  Distance  Comparison  For  The  Three  Guidance  Laws 

(Scenario  #1) 
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2.  Scenario  #2(Varying  Target  Acceleration) 


Figure  5.46  Miss  Distance  Comparison  For  The  Three  Guidance  Laws 

(Scenario  #2) 
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Figure  5.47  Missile  Acceleration  Comparison  For  The  Three  Guidance  Laws 

(Scenario  #2) 
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3.  Discussion 

The  three  dimensional  miss  distance  may  be  improved  by  estimating  the  3-D 
target’s  evasive  maneuvei .  One  way  to  estimate  the  3-D  target  acceleration  is  by  utilizing 
dynamic  image  processing.  Three  setups  were  considered: 

1.  two  cameras, 

2.  a  camera  and  a  radar, 

3.  only  one  camera. 

In  the  single  camera  setup,  actual  range  is  achieved  by  guessing  the  target’s 
physical  dimensions.  The  target’s  3-D  motion  parameters  can  be  estimated  by  utilizing  two 
consecutive  image  frames.  The  target’s  acceleration  may  be  computed  by  taking  the  second 
derivative  after  filtering  the  target’s  motion  data.  Alternatively,  the  target’s  3-D  motion 
may  be  processed  by  Kalman  filters  to  estimate  its  acceleration  components.  This 
additional  information  about  the  target’s  behavior  is  injected  in  suitable  control  laws  to 
improve  the  missile’s  homing  performance. 

Proportional  navigation,  augmented  proportional  navigation,  and  optimal 
guidance  laws  were  derived  for  use  in  a  three  dimensional  environment.  The  classical 
proportional  navigation  guidance  law  tracks  a  target  with  good  accuracy,  especially  if  the 
target  maneuvers  at  long  time  to  go.  However,  when  compared  with  augmented 
PROPNAV  and  optimal  guidance,  PROPNAV  requires  higher  missile  acceleration 
capabilities.  A  plausible  guidance  law  is  one  that  issues  missile’s  commands  proportional 
to  the  miss  distance  that  would  result  if  the  missile  made  no  further  corrections.  Augmented 
proportional  navigation  was  derived  using  this  heuristic  argument.  For  a  constant  target 
maneuver,  augmented  proportional  navigation  increases  the  missile  percentage  of  kill.  For 
a  non  constant  evasive  maneuver,  APROPNAV  does  not  always  guarantees  less  miss 
distance  than  PROPNAV.  However,  APROPNAV  requires  less  missile  acceleration 
capabilities  than  PROPNAV.  Optimal  guidance  was  derived  for  a  missile  with  a  single  lag 
guidance  system.  Optimal  guidance  provides  compensation  for  the  missile's  guidance 
system  dynamics.  The  optimal  guidance  law  requires  the  least  missile  acceleration 
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capability  of  the  three  guidance  laws.  In  fact,  this  law  is  derived  in  order  to  drive  the  miss 
distance  to  zero  while  minimizing  a  performance  index  made  up  of  the  integral  of  the 
square  of  the  control.  Clearly,  the  optimal  guidance  law  presents  the  least  miss  distance  of 
the  three  guidance  laws.  However,  it  requires  a  missile  with  complex  signal  processing 
capabilities.  The  homing  capabilities  of  the  missile  can  be  dramatically  increased  by 
identifying  the  target’s  evasive  maneuver  and  injecting  this  information  into  the 
APROPNAV  (especially  for  a  constant  target  maneuver)  or  optimal  guidance  control 
algorithms.  The  optimal  control  algorithm  guarantees  extraordinary  performance.  Utilizing 
optimal  guidance,  especially  against  highly  responsive  targets,  can  be  the  difference 
between  failure  and  success. 

B.  RECOMEND  ATIONS 

It  is  recommended  to  continue  this  research  by  simulating  the  overall  system  (i.e. 
estimating  the  3-D  target’s  evasive  maneuver  from  two  consecutive  image  frames  and 
injecting  this  data  into  the  tridimensional  missile/target  engagement  simulation  programs 
developed  in  this  thesis).  The  simulations  developed  in  this  thesis  are  very  generic  and 
easily  adapted  to  different  conditions  (i.e.  for  systems  with  different  dynamics  and  initial 
conditions).  The  consequences  of  the  image  measurement  errors  in  die  target  acceleration 
estimation  and  ultimately  in  the  miss  distance  can  be  investigated.  Finally,  it  is 
recommended  that  electronic  counter  measures  (ECCM)  be  added  to  the  target’s  evasion 
maneuver  in  order  to  evaluate  their  effects  on  the  miss  distance. 


APPENDIX  A  -  MISSILE/TARGET  THREE  DIMENSIONAL 
SIMULATION  USING  PROPORTIONAL  NAVIGATION  GUIDANCE 


%  Written  by:  Rui  Manuel  Alves  Francisco 
%  Date:  14  October  1992 

%  This  Program  simulates  the  terminal  phase  of  a  3-D  missile/target 
%  engagement  using  classical  proportional  navigation. 

clear 

clg 

%  DEFINE  CONSTANTS 
dt  =  .01;  %  Sampling  time 
Tf  =  100;  %  maximum  simulation  time 
kmax  =  Tf/dt+1; 
n  =  3;  %  navigation  constant 
N  =  [n  0 
On]; 

%  DEFINE  STATE  EQUATIONS 
%  Missile 

%  Xm  =  [xm  -  raissiie's  x  coordinate 
%  xmd  =  missile's  speed  (x  coordinate) 

%  ym  =  missile's  y  coordinate 

%  ymd  =  missile's  speed  (y  coordinate) 

%  zm  =  missile's  z  coordinate 
%  zmd  =  missile's  speed  (z  coordinate)] 

Am  =  [0  1  0  0  0  0 
000000 
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000  1  00 
000000 
000001 
000000]; 

Bm  =  [000 
100 
000 
0  10 
000 
00  1]; 

%  Target 

%  Xt  =  [xt  =  target’s  x  coordinate 
%  xtd  =  target’s  speed  (x  coordinate) 
%  yt  =  target’s  y  coordinate 
%  ytd  =  target's  speed  (y  coordinate) 
%  zt  =  target's  z  coordinate 
%  ztd  =  target's  speed  (z  coordinate)] 
At  =  [0  1  0  0  0  0 
000000 
000100 
000000 
000001 
000000]; 

Bt  =  [0  0  0 
100 
000 
010 
000 
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0  0  1]; 

%  Seeker  (Rada:-) 

%  Xsk  =  [beta_pitch  =  seeker’s  pitch  angle 
%  betad_pitch  =  seeker’s  pitch  angle  rate 

%  beta_yaw  =  seeker’s  yaw  angle 

%  betad_yaw  =  seeker’s  yaw  angle  rate] 

Ask  =  [  0  1  0  0 

-100  -20  0  0 
0  0  0  1 
0  0  -100-20]; 

Bsk  =  [00 
1000 
0  0 
0  100]; 

%  Guidance  System 

%  Xgs  =  [a_m_pitch  =  missile’s  pitch  acceleration 
%  a_m_yaw  =  missile’s  yaw  acceleration] 

Ags  =  [-10 

0-1]; 

Bgs  =  [10 

01]; 

%  INITIALIZE  STATE  VARIABLES  (when  the  missile  enters  into  the  terminal 
phase  of  flight) 

%  Missile 

Xm(:,l)  =  [  0  %  The  missile  is  in  a  collision  triangle 

2828  %  with  the  target  when  the  missile  enters  into 

0  %  the  terminal  phase  of  flight 

1000 
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0 

47.1339]; 

%  Target 
Xt(:,l)  =  [30000 
0 
0 

1000 

500 

0]; 

%  DISCRETE  REPRESENTATION 
[PHIm.DELm]  =  c2d(Am,Bm,dt); 

[PHIt.DELt]  =  c2d(At,Bt,dt); 

[PHIsk,DELsk]  =  c2d(4sk,Bsk,dt); 

[PHIgs.DELgs]  =  c2d(Ags,Bgs,dt); 

%LINE  OF  SIGHT  (LOS)  INFORMATION.  INITIAL  CONDITIONS. 
%  Missile 

%  LAMBDA_m  =  Missile’s  LOS  from  the  global  coordinate  system 
%  LAMBDA_m  =  [LAMBDA_m_pitch  =  Missile’s  pitch  LOS  angle 
%  LAMB D A_m_yaw  =  Missile’s  yaw  LOS  angle] 

LAMBDA_m(:,l)  =  [atan2(Xni(5,l),sqrt(Xm(Ll)A2+Xm(3,l)A2)); 
atan2(Xm(3,l),Xm(l,l))]; 

%  Target 

%  LAMBDA_t  =  Target’s  LOS  from  the  global  coordinate  system 
%  LAMBDA_t  =  [LAMBPA_t_pitch  =  Target’s  pitch  LOS  angle 
%  LAMB D A_t_y aw  =  Target’s  yaw  LOS  angle] 

LAMBDA_t(:,l)  =  [atan2(Xt(5,l),sqrt(Xt(U)A2+Xt(3,l)A2)); 

atan2(Xt(3,l),Xt(l,l))]; 

%  LOS  from  Missile  to  Target 

1(M 


%  LAMBDA  =  LOS  from  Missile  to  Target. 

%  LAMBDA  =  [LAMBDA_pitch  =  LOS  angle  in  pitch 
%  LAMBDA_yaw  =  LOS  angle  in  yaw], 

LAMBDA(:,1)  =  [atan2((Xt(51l)-Xm(5,l)),sqrt((Xt(l,l)-Xm(l,l))A2 
+(Xt(3,l)-Xm(3,l))A2)); 

atan2((Xt(3,l)-Xm(3,l)),abs(Xt(l,l)-Xm(l,l)))]; 

%  MISSILE  and  TARGET  FLIGHT  PATH  ANGLES  INFORMATION 
%  Missile 

%GAMMA_m  =  [GAMMA_m_pitch  =  Missile’s  flight  path  angle  in  pitch 
%  GAMMA_m_yaw  =  Missile’s  flight  path  angle  in  yaw] 

GAMMA_m(:,l)  =  [atan2(Xm(6,l),sqrt(Xm(2,l)A2+Xm(4,l)A2)); 
atan2(Xm(4,l),Xm(2,l))]; 

%  Target 

%GAMMA_t  =  [GAMMA_t _pitch  =  Target’s  flight  path  angle  in  pitch 
%  GAMMA_t_YAW  =  Target’s  flight  path  angle  in  yaw] 

GAMMA_t(:,l)  =  [atan2(Xt(6,l),sqrt(Xt(2,l)A2+Xt(4,l)A2)); 

atan2(Xt(4, 1 )  ,Xt(2, 1 ))] ; 

%  RANGE  INFORMATION 
%  Missile 

%  Rm  =  Missile’s  range 

Rm(l)  =  sqrt(Xm(l,l)A2  +  Xm(3,l)A2  +  Xm(5,l)A2); 

%  Target 

%  Rt  =  Target’s  range 

Rt(l)  =  sqrt(Xt(l,l)A2  +  Xt(3,l)A2  +  Xt(5,l)A2); 

%  Missile/Target  relative  distance 
%  R  =  [Rmtx  =  Missile/Target  x  coordinate  range 
%  Rmty  =  Missile/Target  y  coordinate  range 

%  Rmtz  =  Missile/Target  z  coordinate  range 
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%  Rmt  =  Missile/Target  relative  distance(Rmt=sqrt(RmtxA2+RmtyA2+RmtzA2))] 
R(:,l)  =  [Xt(l,l)-Xm(l,l) 

Xt(3,l)-Xxn(3,l) 

Xt(5,l)-Xm(5,l) 

sqrt((Xt(  1 , 1  )-Xm(  1 , 1  ))A2+(Xt(3, 1  )-Xm(3, 1  ))A2+(Xt(5, 1  )-Xm(5, 1 ))  A2)] ; 

%  VELOCITY  INFORMATION 
%  Missile 

%  Vm  =  Missile’s  Speed 

Vm(l)  =  sqrt(Xm(2, 1)A2+Xm(4, 1)A2+Xm(6, 1)A2); 

%  Target 

%  Vt  =  Target's  Speed 

Vt(l)  =  sqrt(Xt(2, 1  )A2+Xt(4, 1  )A2+Xt(6, 1  )A2); 

%  Speed  along  the  pitch  and  yaw  LOS.  Pitch  and  yaw  closing  speeds 
VLpitch(l)  =  Vt(l)*cos(LAMBDA(2,l)-GAMMA_t(2,l)); 

Vm_pitch(l)  =  Vm(l):t‘cos(LAMBDA(2)l)-GAMMA_m(2,l)); 

Vc_pitch(l)  =  Vm_pitch(l)*cos(GAMMA_m(  1 ,1)-LAMBDA(  1,1)) 

-  Vt.4?itch(l)>,,cos(GAMMA_t(  1,1)-LAMBDA(  1,1)); 

Vt_yaw(l)  =  Vt(  l)*cos(GAMMA_t(  1,1)); 

Vm_yaw(l)  =  Vm(l)*cos(GAMMA_m(l,l)); 

Vc_yaw(l)  =  Vm_yaw(l)’!<cos(GAMMA_m(2,l)-LAMBDA(2,l)) 
-Vt_yaw(l),|‘cos(GAMMA_t(2,l)~LAMBDA(2,l)); 

Vc  =  [Vc_pitch(l)  0 

0  Vc_yaw(l)3, 

%  SEEKER  and  GUIDANCE  SYSTEM  INITIAL  CONDITIONS  and  INPUTS. 
%  Seeker 
Xsk(:,l)  =  [0 
0 
0 
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0]; 


U$k(:,l)  =  LAMBDA(:,1);  %  Seeker  input 
%  Guidance  System 
Xgs(:,l)  =  [0 
0]; 

Ugs(:,l)  =  N*Vc*[Xsk(2,l)  %  Guidance  system  input 
Xsk(4,l)]; 

%  TIME 

%  Time  =  Time  vector 
TIME(l)  =  0; 

%  Tgo  =  Time  to  go 
Tgo(l)  =  R(4,l)/Vc_piich(l); 

%  SIMULATE  THE  SYSTEM 
for  ti  =  0:. 25:2 1.25 
for  i=  l:kmax-l 

%  Calculate  components  of  the  missile’s  pitch  acceleration 
a_mx_pitch(i)  =  -(Xgs(U)*sin(LAMBDA(U))*cos(LAMBDA(2,i))); 
a_my_pitch(i)  =  -(Xgs(l,i)*sin(LAMBDA(l,i))’‘sin(LAMBDA(2,i))); 
a_mz_pitch(i)  =  Xgs(l,i)*cos(LAMBDA(l,i)); 

%  Calculate  components  of  the  missile’s  yaw  acceleration 
a_mx_yaw(i)  =  -Xgs(2,i)*$in(LAMBDA(2,i)); 
a_my_yaw(i)  =  Xgs(2,i)*cos(LAMBDA(2,i)); 

%  Compute  overall  missile  acceleration 
a_mx(i)  =  a_mx_pitch(i)+a_mx_yaw(i); 
a_my(i)  =  a_my_pitch(i)+a_my_.yaw(i); 
a_mz(i)  =  a_mz_pitch(i); 
am(:,i)  =  [a_mx(i) 
a_my(i) 
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a_mz(i)]; 


%  Compute  missile’s  acceleration  magnitude 
a_m(i)  =  sqrt(a_mx(i)A2  +  a_my(i)A?  +  a_mz(i)A2); 

%  Generate  target’s  evasive  maneuver  (we  assume  that  these  accelerations,  along 
%  the  three  cartesian  axis,  are  estimated  using  the  missile’s  image  processing 
%  capabilities) 

if  TIME(i)  >=  ti/2  %  target  starts  evasive  maneuver 
a_tx(i)  =  3*32,2*sin(GAMMA_t(2,i)); 
a_ty(i)  =  4*32,2*cos(GAMMA_t(2,i)); 
a_tz(i)  =  3*32.2*cos(GAMMA_t(l,i)); 
else 

a_tx(i)  =  0.0; 
a_ty(i)  =0.0; 
a_tz(i)  =  0.0; 
end 

at(:,i)  *  |a_tx(i) 
a_ty(i) 
a  _tz(i)]; 

%  Compute  magnitude  of  the  target’s  acceleration 
a_t(i)  =  sqrt(a_tx(i)A2  +  a_ty(i)A2  +  a_tz(i)A2); 

%  Update  missile  states 

Xm(:,i+1) »  PHIm*Xm(:,i)+DELm*am(:,i); 

%  Update  target  states 

Xt(:,i+1)  =  PHIt*Xt(:,i)+DELt*at(:,i); 

%  Update  seeker  states 

Xsk(:,i+1)  =  PHIsk*Xsk(:,i)+DELsk*Usk(:,i); 

%  Update  Guidance  System  states 
Xgs(:,i+1)  =  PHIgs*Xgs(:,i)+DELgs*Ugs(:,i); 
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%  Limit  yaw  and  pitch  accelerations  to  25  g’s 
if  ab$(Xgs(l,i+l))  >  805,0 
Xgs(l,i+1)  =  805.0  *sign(Xgs(t,i+l)); 
end 

if  abs(Xg$(2,i+l))  >  805.0 
Xgs(2,i+1)  =  805.0  *sign(Xgs(2,i+l)); 
end 

%  Update  LOS  angles 

LAMBDA(:,i+l)  =  [atan2((Xt(5,i+i)-Xm(5,i+l)),sqrt((Xm(l,i+l)-Xt(Li+l))A2 
+(Xm(3,i+l)-Xt(3,i+l))A2)); 

atan2((Xt(3,i+l)-Xm(3,i+l)),(abs(Xm(l,i+l)-Xt(l,i+l))))]; 

Usk(:,i+1)  =  LAMBDA(:,i+l); 

LAMBDA_m(:,i+l)  =  [atan2(Xm(5,i+l),sqrt(Xm(l,i+-l)A2+Xm(3fi+l)A2)); 

atan2(Xm(3,i+l  ),Xm{  l,i+ 1 ))] ; 

LAMBDA_t( :  ,i+ 1 )  =  [atan2(Xt(5,i+ 1  ),sqrt(Xt(  1  ,i+ 1  )A2+Xt(3,i+ 1 )  A2)) ; 

atan2(Xt(3,i+ 1  ),Xt(l  ,i+ 1 ))] ; 

%  Update  flight  path  angles 

GAMMA_m(:,i+l)  =  [atan2(Xm(6,i+l),sqrt(Xm(2.i+l)A2-fXm(4,i+l)A2)); 

atan2(Xm(4,i+ 1  ),Xm(2,i+ 1 ))]; 

GAMMA_t(:,i+l)  =  [atan2(Xt(6,i+l),sqrt(Xt(2,i+l)A2+Xt(4,i+l)A2)); 

atan2(Xt(4,i+l),Xt(2,i+l))]; 

%  Update  Range  Information 

Rm(i+1)  =  sart(Xm(l,i+l)A2  +  Xm(3,i+1)A2  +  Xm(5a+l)A2); 

Rt(i+1)  =  sqrt(Xt(l,i+l)A2  +  Xt(3,i+1)A2  +  Xt(5,i+1)A2); 

R(:,i+1)  =  [Xt(l,i+1)-Xm(3,i+1); 

Xt(3,i+1)-Xm(3,i+1); 

Xt(5,i+l)-Xm(5,i+i); 

sqrt((Xt(  1  ,i+ 1  )-Xm(  1  ,i+ 1  ))A2+(Xt(3,i+ 1  )-Xm(3,H- 1 ))  A2 
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+(Xt(5,i+l)-Xm(5,i+l))A2)]; 


%  Update  Velocity  Information 

Vm(i+1)  =  sqrt(Xm(2,i+ 1)  A2+Xm  (4,i+ 1  )A2+Xm(6,i+  1)A2); 

Vt(i+1)  =  sqrt(Xt(2,i+ 1  )A2+Xt(4,i+1  )A2+Xt(6,i+ 1  )A2); 

Vt_pitch(i+1)  =  Vt(i+l)*cos(LAMBDA(2,i+ 1  )-GAMMA_t(2,i+l )); 
Vm_pitch(i+1)  =  Vm(i+l)*cos(LAMBDA(2,i+l)-GAMMA_m(2,i+l)); 
Vc_pitch(i+1)  =  Vm  j3itch(H4)*cos(GAMMA_m(l,i+l)-LAMBDA(l,i+l)) 
-  V  t_pitch(i+ 1  )*cos(GAMMA_t(  1  ,i+  1)-LAMBD  A(  1,1+1)); 
Vt_yaw(i+1)  =  Vt(i+ 1  )*cos(G  AMMA_t(  1  ,i+ 1 )); 

Vm_yaw(i+1)  =  Vm(i+l)*cos(GAMMA_m(l,i+l)); 

Vc_yaw(i+1)  =  Vm_yaw(i+l)*cos(GAMMA_m(2,i+l)-LAMBDA(2,i+l)) 
-Vt_yaw(i+l)*cos(GAMMA_t(2,i+l)-LAMBDA(2,i+l)); 
Vc  =  [Vc_pitch(i+1)  0 

0  Vc_yaw(i+1)]; 

%  Update  guidance  system  input 
Ugs(:,i+1) «  N*Vc*[Xsk(2,i+l) 

Xsk(4,i+1)]; 

%  Update  Time/time  to  go 
TTME(i+l)  =TIME(i)+dt; 

Tgo(i+l)  =  R(4,i+ 1  )/Vc_pitch(i+ 1 ); 

%  Check  for  closest  point 
if  (R(4,i)  <  R(4,i+1)), break, end 
end; 

%  Save  information  for  plotting  and  evaluation 
Rl(4*ti+1)  =  R(4,i);  %  miss  distance 
Ti(4*ti+1)  =  ti/2;%  starting  time  of  evasive  maneuver  (EM) 
tgo(4*ti+l)  =  i*dt-ti/2;  %  time  to  go  until  end  of  flight 
if  ti  ==  12.0  %  Record  information  for  a  target  that 
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%  initialized  the  evasive  maneuver  6  sec  after  the  missile 
%  entered  into  the  terminal  phase  of  flight. 

TGO  =  tgo(49); 

Xseeker  =  Xsk(:,l:i); 

Xgsys  =  Xgs(:,l:i); 
lambda_m  =  LAMBDA_m(:,l:i); 
lambda_t  =  LAMBDA_t(:,l:i); 
lambda  =  LAMBDA(:,l:i); 
gamma_rn  =  GAMMA_m(:,l:i); 
gamma_t  =  GAMMA_t(:,l:i); 
r  =  R(:,l:i); 
vm  =  Vm(l:i); 
vt  =  Vt(l:i); 

vm_pitch  =  Vm_pitch(l:i); 
vt_pitch  =  Vt_pitch(l:i); 
vm_yaw  =  V  m_yaw(  l:i); 
vt__yaw  =  Vt_yaw(l:i); 
vc_pitch  =  Vc_pitch(l:i); 
vc_yaw  =  Vc_yaw(l:i); 
tGO  =Tgo(l:i); 
a_M  =  am(:,l:i); 
a_T  =  at(:,l:i); 

A_t =  a_t(l:x); 

A_m  =  a_m(l:i); 
time  =  TIME(l:i); 
end 

clear  R; 

R(:,l)  =  [Xt(l,l)-Xm(l,l); 
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Xt(3,l)-Xm(3,l); 

Xt(5,l)-Xm(5,l); 

sqrt((Xt(l,l)-Xm(l,l))A2+(Xt(3,l)-Xm(3,l))A2+(Xt(5,l)-Xm(5,l))A2)]; 

end; 

save  thesis  lp343  R1  tgo  Ti  tGO  missile  target  TGO  Xseeker  Xgsys  lambda_m 
lambda_t  lambda  gamma_m  gamma_t  r  vm  vt  vm_pitch  vm_yaw  vt_pitch 
vt_yaw  vc_pitch  vc_yaw  a_M  a_T  time  A_t  A_m 
PLOTS 

%  Miss  distance  information 

plot(Ti,R  1  ),title( ‘ MISS  DISTANCE  vs  INITIAL  TIME,  PROPNAV’) 

xlabelO INITIAL  TIME  -  SEC’),ylabel(‘MISS  -  FEET’) 

print  -dps  Rlapl 

Ipstoepsi  Rlapl.ps  Rlapl.epsi 

pause,clg 

plot(tgo,R  l),title(‘MISS  DISTANCE  vs  TIME  TO  GO,  PROPNAV’) 

xlabelCTEME  TO  GO  -  SEC).ylabel(‘MISS  -  FEET’) 

print  -dps  Rlbpl 

Ipstoepsi  Rlbpl.ps  Rlbpl.epsi 

pause, clg 

%  Missile  acceleration  information 

plot(time,A_m),title(‘MISSILE  ACCELERATION  MAGNITUDE  vs  TIME, 

PROPNAV’) 

Alauw\  i±±VklJi-v Ij\_  ruDi/oiiv  x  ) 

print  -dps  A__mpl 

Ipstoepsi  A_mpl.ps  A_mpl.epsi 

pause,clg 

plot(time,Xgsys(l,:)),title(‘MISSILE  PITCH  ACCELERATION  vs  IIME, 

PROPNAV’) 
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xlabel(‘TTME  -  SF.C’),ylabei(‘FEET/SECA2’) 
print  -dps  Xgsyslpl 
Ipstoepsi Xgsyslpl.ps  Xgsyslpl.epsi 
pause, clg 

piot(time,Xgsys(2,:)),title(‘MISSILE  YAW  ACCELERATION  vs  TIME, 

PROPNAV’) 

xiabel(‘TTME  -  S  EC  ’ )  ,y  label(  ‘  FEET/S  EC  A2  ’ ) 
print  -dps  Xgsys2pl 
Ipstoepsi Xgsys2pl.ps  Xgsys2pl.epsi 
pause, clg 

%  Target  acceleration  information 

p lot( time, A_t),title( ‘TARGET  ACCELERATION  MAGNITUDE  vs  TIME, 
PROPNAV’) 

xlabelHTME  -  SEC’),ylabel(‘FEET/SECA2’) 

print  -dps  A_tpl 

Ipstoepsi  A_tp  1  .ps  A_tp  1  .epsi 

pause.clg 

%  Seeker  pitch  and  yaw  angles 

plot(time,Xseeker(l,:)),titie(‘SEEKER  PITCH  ANGLE  vs  TIME,  PROPNAV’) 

xlabel(‘TIME  -  SEC’),ylabel(‘RAD’) 

print  -dps  Xseekerlpl 

Ipstoepsi  Xseekerlpl.ps  Xseekerlpl. epsi 

pause.clg 

plot(time,Xseeker(3,:)),title(‘SEEKER  YAW  ANGLE  vs  TIME,  PROPNAV’) 

xlabel(‘TIME  -  SEC’),ylabel(‘RAD’) 

print  -dps  Xseeker2pi 

Ipstoepsi  Xseeker2pl.ps  Xseeker2pl.epsi 

pause.clg 
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plot(time,Xseeker(2,:)),title(‘SEEKER  PITCH  ANGLE  RATE  vs  TIME, 

PROPNAV’) 

xiabei(‘TIME  -  SEC’),ylabel(‘RAD/SEC*) 
print  -dps  Xseeker3pl 
Ipstoepsi  Xseeker3pl.ps  Xseeker3p  1  .epsi 
pause.clg 

plot(time,Xseeker(4,:)) , title(  ‘  S  EEKER  YAW  ANGLE  RA'IE  vs  TIME, 

PROPNAV’) 

xlabel(‘TIME  -  SEC’),yiabel(‘RAD/SEC’) 
print  -dps  Xseeker4pl 
Ipstoepsi  Xseeker4pl.ps  Xseeker4pl.epsi 
pause.clg 

%  Range  information 

plot(time/(4,:)),title(‘RANGE  vs  TIME,  PROPNAV’) 

xlabel(‘TIME  -  SEC’),yiabel(‘FEET’) 

print  -dps  rpl 

Ipstoepsi  rpl.ps  rpl. epsi 

pause.clg 

%  Missile  velocity  information 

plot(time,vm),iitle(‘MISSILE  SPEED  vs  TIME,  PROPNAV’) 

xlabel(‘TIME  -  SEC’),ylabel(‘FEET/SEC’) 

print,  -dps  vmpl 

Ipstoepsi  vmpl.ps  vmpl. epsi 

pause.clg 

%  Target  velocity  information 

plot(time,vt),title( ‘TARGET  SPEED  vs  TIME,  PROPNAV’) 
xlabel(‘TIME  -  SEC’),ylabel(‘FEET/SEC’) 
print  -dps  vtpl 
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Ipstoepsi  vtpl.ps  vtpl.epsi 
pause, clg 

%  Closing  velocity  information 

plot(time,vc_pitch),title( ‘PITCH  CLOSING  SPEED,  PROPNAV’) 

xlabel(‘TIME  -  SEC’).ylabel(‘FEET/SEC’) 

print  -dps  vclpl 

Ipstoepsi  vclpl.ps  vclpl.epsi 

pause, clg 

plot( time, vc_yaw),title(‘ YAW  CLOSING  SPEED  vs  TIME,  PROPNAV’) 

xlabelCTIME  -  SECXylabelCFEEl/SEC’) 

print  -dps  vc2pl 

Ipstoepsi  vc2pl.ps  vcipl.epsi 

pause.clg 


APPENDIX  B  -  MISSILE/TARGET  THREE  DIMENSIONAL 
SIMULATION  USING  AUGMENTED  PROPORTIONAL  NAVIGATION 
GUIDANCE 


%  Written  by:  Rui  Manuel  Alves  Francisco 
%  Date:  10  November  1992 

%  This  Program  simulates  the  terminal  phase  of  a  3-D  missile/target 
%  engagement  using  augmented  proportional  navigation  guidance. 

clear 

clg 

%  DEFINE  CONSTANTS 
dt  =  .01;  %  Sampling  time 
Tf  =  100;  %  maximum  simulation  time 
km/  r  =  Tf/dt+1; 
n  =  3;  %  navigation  constant 
N  =  [n  0 
On]; 

%  DEFINE  STATE  EQUATIONS 
%  Missile 

%  Xm  =  [xm  -  missile's  x  coordinate 
%  xmd  =  missile's  speed  (x  coordinate) 

%  ym  =  missile's  y  coordinate 

%  ymd  =  missile's  speed  (y  coordinate) 

%  zm  =  missile's  z  coordinate 

%  zmd  =  missile's  speed  (z  coordinate)] 

Am  =  [0  1  0  0  0  0 
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000000 
000 100 
000000 
00000  1 
0  0  0  0  0  0], 

Bm  =  [0  0  0 
100 
000 
0  10 
000 
0  0  1]; 

%  Target 

%  Xt  =  [xt  =  target’s  x  coordinate 
%  xtd  =  target’s  speed  (x  coordinate) 
%  yt  =  target’s  y  coordinate 
%  ytd  =  target's  speed  (y  coordinate) 
%  zt  =  target's  z  coordinate 
%  ztd  =  target's  speed  (z  coordinate)] 
At  =  [0  1  0  0  0  0 
000000 
000100 
000000 
000001 
000000]; 

Bt  =  [0  0  0 
100 
000 
0  1  0 


117 


0  0  1]; 

%  Seeker  (Radar) 

%  Xsk  =  [beta .pitch  =  seeker’s  pitch  angle 
%  betad_pitch  =  seeker’s  pitch  angle  rate 

%  beta.yaw  =  seeker’s  yaw  angle 

%  betad.yaw  =  seeker’s  yaw  angle  rate] 

Ask  =  [  0  1  0  0 

-100  -20  0  0 

0  0  0  1 
0  0  -100-20]; 

Bsk  =  [00 
100  0 
0  0 
0  100]; 

%  Guidance  System 

%  Xgs  =  [a_m_pitch  =  missile’s  pitch  acceleration 
%  a_m_yaw  =  missile’s  yaw  acceleration] 

Ags  =  [-1  0 


0-1]; 


Bgs  =  [10 


%  INITIALIZE  STATE  VARIABLES  (when  the  missile  enters  into  the  terminal 
phase  of  flight) 


%  Missile 
Xm(:,l)  =  [  0 


%  The  missile  is  in  a  collision  triangle 
%  with  the  target  when  the  missile  enters  into 
%  the  terminal  phase  of  flight 
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Xt(:.l)  =  [30000 
0 
0 

1000 

500 

0]; 

%  DISCRETE  REPRESENTATION 
[PHIm.DELm]  =  c2d(Am,Bm,dt); 

[PHHt,DELt]  =  c2d(At,Bt,dt); 

[PHIsk.DELsk]  =  c2d(Ask,Bsk,at); 

[PHIgs.DELgs]  =  c2d(Ags,Bgs,dt); 

%  LINE  OF  SIGHT  (LOS)  INFORMATION.  INITIAL  CONDITIONS. 
%  Missile 

%  LAMBDA_m  =  Missile’s  LOS  from  the  global  coordinate  system 
%  LAMBDA_m  =  [LAMBDA_m_pitch  =  Missile’s  pitch  LOS  angle 
%  LAMB D A_m_y aw  =  Missile’s  yaw  LOS  angle] 

LAMBDA_m(:,l)  =  [atan2(Xm(5,l),sqrt(Xm(l,l)A2+Xm(3,l)A2)); 
atan2(Xm(3,l),Xm(l,l))]; 

%  Target 

%  LAMEDA_t  =  Target’s  LOS  from  the  global  coordinate  system 
%  LAMBDA_t  =  [LAMBDA_t_pitch  =  Target’s  pitch  LOS  angle 
%  LAMB  D  A_t_y aw  =  Tai'get’s  yaw  LOS  angle] 

LAMBDA_t(:,l)  =  [atan2(Xt(5,l),sqrt(Xt(l,l)A2+Xt(3,l)A2)); 
atan2(Xt(3, 1  ),Xt(  1,1))]; 


%  LOS  from  Missile  to  Target 
%  LAMBDA  =  LOS  from  Missile  to  Target. 

%  LAMBDA  =  [LAMBDA_pitch  =  LOS  angle  in  pitch 
%  LAMBDA_yaw  =  LOS  angle  in  yaw]; 

LAMBDA(:,1)  =  [atan2((Xt(5,l)-Xm(5.1)),sqrt((Xt(l,l)-Xm(U))A2 
+(Xt(3,l)-Xm(3,l))A2)); 

atan2((Xt(3,l)-Xm(3,l)),abs(Xt(l,l)-Xm(l,i)))]; 

%  MISSILE  and  TARGET  FLIGHT  PATH  ANGLES  INFORMATION 
%  Missile 

%GAMMA_m  =  [GAMMA_m_pitch  =  Missile’s  flight  path  angle  in  pitch 
%  GAMMA_m_yaw  =  Missile’s  flight  path  angle  in  yaw] 

GAMMA_m(:,l)  =  [atan2(Xm(6,l).sqrt(Xm(2,l)A2+Xm(4.1)A2)); 
atan2(Xm(4, 1 )  ,Xm(2, 1 ))] ; 

%  Target 

%GAMMA_t  =  [GAMMA_t_pitch  =  Target’s  flight  path  angle  in  pitch 
%  G AMMA_t_ Y AW  -  Target’s  flight  path  angle  in  yaw] 

GAMMA_t(:,l)  =  [atan2(Xt(6,l),sqrt(Xt(2(l)A2+Xt(4,l)A2)); 

atan2(Xt(4,  l  ),Xt(2, 1 ))] ; 

%  RANGE  INFORMATION 
%  Missile 

%  Rm  -  Missile’s  range 

Rm(l)  «  sqrt(Xm(l,l)A2  +  Xm(3,l)A2  +  Xm(5,i)A2); 

%  Target 

%  Rt  =  Target’s  range 

Rt(l)  =  sqrt(Xt(l,l)A2  +  Xt(3,l)A2  +  Xt(5,l)A2); 

%  Missile/Target  relative  distance 
%  R  =  [Rmtx  =  Missile/Target  x  coordinate  range 
%  Rmty  =  Missile/Target  y  coordinate  range 
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%  Rmtz  =  MissileA'arger  z  coordinate  range 

%  Rmt  =  MissileA'arget  relative  distance(Rmt=sqrt(RrntxA2+RmtyA2+RmtzA2))] 
R(:,l)  =  [Xt(l,l)-Xm(l,l) 

Xt(3,l)-Xm(3.1) 

Xt(5,l)-Xm(5,l) 

sqrt((Xt(l,l)-Xm(U))A2+(Xt(3,l)-Xm(3fl))A2+(Xt(5.1)-Xm(5,l))A2)]; 

%  VELOCITY  INFORMATION 
%  Missile 

%  Vm  =  Missile’s  Speed 

Vm(l)  =  sqrt(Xm(2,l)A2+Xm(4,l)A2+Xin(6,l)A2); 

%  Target 

%  Vt  -  Target’s  Speed 

Vt(l)  =  sqrt(Xt(2, 1  )A2+Xt(4, 1  )A2+Xt(6,  i  )A2); 

%  Speed  along  the  pitch  and  yaw  LOS.  Pitch  and  yaw  closing  speeds 
Vt_pitch(l)  =  Vt(l)*cos(LAMBDA(2ll)-GAMMA_t(2,l)); 

Vm_pitch(l)  =  Vm(l)*cos(LAMBDA(2,l)-GAMMA_m(2,l)); 

Vc_pitch(l)  =  Vm_pitch(l)*cos(GAMMA_m(l,l)-LAMBDA(l,l)) 

-  Vt_pitch(  1  )*cos(G  AMMA_t(  1 , 1)-LAMB  D  A(  1 , 1 )) ; 

Vt_yaw(l)  =  Vt(T)*cos(GAMMA_t(l,l)); 

Vm_yaw(l)  =  Vm(l)*cos(GAMMA_m(l,l)); 

Vc_yaw(l)  =  Vm_yaw(l)*cos(GAMMA_m(2,l)-LAMBDA(2,l)) 
-Vt_yaw(l)'l<cos(GAMMA_t(2,l.)-LAMBDA(2,l)); 

Vc  =  [Vc_pitch(l)  0 

0  Vc_yaw(l)]; 

%  SEEKER  and  GUIDANCE  SYSTEM  INITIAL  CONDITIONS  and  INPUTS. 
%  Seeker 
Xsk(:,l)  =  [0 
0 
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IS 


tn 


as 


0 

0]; 

Usk(:,l)  =  LAMBDA(:,1);  %  Seeker  input 
%  Guidance  System 
Xgs(:.l)  =  [0 
0]; 

Ugs(:,l)  =  N*Vc*[Xsk(2,l)  %  Guidance  system  input 
Xsk(4,l)]; 

%  Initial  conditions  of  the  target’s  acceleration 

ajtx(l)  =  0; 

a_ty(l)  =0; 

a_tz(l)  =0; 

a_t(l)  =  0; 

a_t_  pitch(l)  =0; 

a_t_yaw(l)  =0; 

TETA_t(  i)  =  0;  %  angle  between  the  acceleration  vector  and  the  yaw  plane 
PHI_t(l)  =  0;  %  yaw  angle  of  the  target’s  acceleration  vector 

%  TIME 

%  Time  =  Time  vector 
TIME(l)  =  0; 

%  Tgo  =  Time  to  go 
Tgo(l)  =  R(4,l)/Vc_pitch(l); 

%  SIMULATE  THE  SYSTEM 
for  ti  =  0:. 25:2 1.25 
for  i  =  l:kmax-l 

%  Calculate  components  of  the  missile’s  pitch  acceleration  vector 
a_mx  pitch(i)  =  -(Xgs(l,i)*sin(LAMBDA(l,i))*cos(LAMBDA(2,i))); 
a_my_pitch(i)  =  -(Xgs(l,i)*sin(LAMBDA(l,i))*sin(LAMBDA(2,i))); 
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a_mz_pitch(i)  =  Xgs(l,i)*cos(LAMBDA(l,i)); 

%  Calculate  components  of  the  missile’s  yaw  acceleration  vector 
a_mx_yaw(i)  =  -Xgs(2,i)*sin(LAMBDA(2,i)); 
a_my_yaw(i)  =  Xgs(2,i)*cos(LAMBDA(2,i)); 

%  Compute  overall  missile  acceleration  components 
a_mx(i)  =  a_mx_pitch(i)+a_mx_yaw(i); 
a_my(i)  =  a_my_pitch(i)+a_my_yaw(i); 
a_mz(i)  -  a_mz_pitch(i); 
am(:,i)  =  [a_mx(i) 
a_my(i) 
a_mz(i)]; 

%  target  acceleration  vector 
at(:,i)  =  [a_tx(i) 
a„ty(i) 
a_tz(i)]; 

%  Compute  magnitude  of  the  missile’s  acceleration 
a_m(i)  =  sqrt(a„mx(i)A2  4-  a_my(i)A2  +  a_mz(i)A2); 

%  Generate  target’s  evasive  maneuver  (we  assume  that  these  accelerations,  along 
%  the  three  cartesian  axis,  are  estimated  using  the  missile’s  image  processing 
%  capabilities) 

if  TIME(i)  >=  ti/2  %  target  starts  evasive  maneuver 
a_tx(i+l)  =  3*32.2*sin(GAMMA_t(2,i)); 
a_ty(i+l)  o  4*32.2*cos(GAMMA_t(2,i)); 
a_tz(i+l)  =  3*32.2,|<cos(GAMMA_t(l,i)); 
else 

a_tx(i+l)  =  0.0; 
a_ty(i+l)  =  0.0; 
a_tz(i+l)  =  0.0; 
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end 

%  Compute  magnitude  of  the  target’s  acceleration 
a_t(i+l)  =  sqrt(a_tx(i+l)A2  +  a_ty(i+l)A2  +  a_tz(i+l)A2); 

%  Update  missile  states 

Xm(:.i+1)  =  PHIm*Xm(:  j)+DELm*am(:,i); 

%  Update  target  states 

Xt(:,i+1)  =  PHIt*Xt(:,i)+DELt*at(:,i); 

%  Update  flight  path  angles 

GAMMA_m(:,i+l)  =  [atan2(Xm(6,i+l),sqrt(Xm(2,i+l)A2+Xm(4,i+l)A2)); 

atan2(Xm(4,i+  l),Xm(2,i+ 1))] ; 

GAMMA_t(:,i+l)  =  [atan2(Xt(6,i+i),sqrt(Xt(2,i+l)A2+Xt(4.i+l)A2)); 

atan2(Xt(4.i+ 1  ),Xt(2,i+ 1 ))] ; 

%  Update  seeker  states 

Xsk(:.i+1)  =  PHIsk*Xsk( :  ,i)+DELsk*Usk( :  ,i); 

%  Update  Guidance  System  states 
Xgs(:,i+1)  =  PHIgs*Xgs(;d)+DELgs*Ugs(:,i); 

%  Limit  yaw  and  pitch  accelerations  to  25  g’s 

if  abs(Xgs(l,i+l))  >  805.0 
Xgs(l,i+1)  =  805.0  *sign(Xgs(l,i+l)); 
end 

if  abs(Xgs(2,i+l))  >  805.0 
Xgs(2,i+1)  =  8G5.0  *sign(Xgs(2,i+l)); 
end 

%  Update  LOS  angles 

LAMBDA(:,i+l)  =  [atan2((Xt(5.i+ 1  )-Xm(5,i+l)),sqrt((Xm(  1  ,i+ 1  )-Xt(l  ,i+ 1  ))A2 
+(Xm(3,i+l)-Xt(3,i+l))A2)); 

atan2((Xt(3,i+l  )-Xm(3,i+  l)),(abs(Xm(  1  ,i+ 1  )-Xt(  1  ,i+l ))))]; 
Usk(:,i+1)  =  LAMBDA(:,i+l); 
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LAMBDA_m(:,i+l)  =  [atan2(Xm(5,i+l),sqrt(Xm(l,i+l)A2+Xm(3,i+l)A2)); 

atan2(Xm(3,i+ 1  ).Xm(  1  ,i+ 1 ))] ; 

LAMBDA_t(:  ,i+ 1)  =  [atan2(Xt(5  ,i+ 1  ).sqrt(Xt(  l,i+ 1  )A2+Xt(3,i+ 1  )A2)); 

atan2(Xt(3  ,i+ 1  ),Xt(  1  ,i+ 1 ))] ; 

%  Update  Range  Information 

Rm(i+1)  =  sqrt(Xm(l,i+l)A2  +  Xm(3,i+1)A2  +  Xm(5,i+1)A2); 

Rt(i+1)  =  sqrt(Xt(l,i+l)A2  +  Xt(3,i+1)A2  +  Xt(5,i+1)A2); 

R(:,i+1)  =  [Xt(l,i+1)-Xm(l,i+1); 

Xt(3  ,i+ 1  )-Xm  (3  ,i+ 1 ) ; 

Xt(5  ,i+  i  )-Xm(5  ,i+ 1 ), 

sqrt((Xt(  l,i+ 1 )  •  Xxn(  1  ,i+l  ))A2+(Xt(3,i+ 1  )-Xm(3.i+ 1 ))  A2 
+(Xt(5,i+lj-Xm(5,i+l))A2)]; 

%  Update  Velocity  Information 

Vm(i+1)  -  sqrt(Xm(2,i+l)A2+Xrn(4,i+l)A2+Xm(6,i+l)A2); 

Vt(i+1)  =  sqrt(Xt(2,i+ 1  )A2+Xt(4,i+ 1  )A2+Xt(6,i+ 1 )  A2) ; 

Vt_pitch(i+1)  =  V  t(i+ 1  )*cos(LAMBDA(2,i+  l)-GAMMA_t(2,i+ 1 )); 
Vm_pitch(i+1)  =  Vm(i+l)*cos(LAMBDA(2,i+l)-GAMMA_m(2,i+l)); 
Vc_pitch(i+1)  =  Vm_p.itch(i+l)*cos(GAMMA_m(i,i+l)-LAMBDA(l,i+l)) 
-Vt_pitch(i+l)*cos(GAMMA_t(  1  ,i+  D-LAMBD  A(  I  ,i+ 1 )); 
Vt_yaw(i+1)  =  Vt(i+l)*cos(GAMMA_t(U+l)); 

Vm_yaw(i+1)  =  Vm(i+l)*cos(GAMMA_m(l,Hl)); 

Vc_yaw(i+1)  =  Vm_yaw(i+l)*cos(GAMMA_m(2,i+l)'LAMBDA(2,i+l)) 

-  Vt_yaw(i+ 1  )*cos(GAMMA_t(2,i+ 1  )-LAMBD  A(2,i+ 1 )); 
Vc  =  [Vc_pitch(i+1)  0 

0  Vc_yaw(i+1)]; 

%  Calculate  angles  of  the  target’s  acceleration 
TETA_t(i+l)  =  atan2(a_tz(i+l),sqrt(a_tx(i+l)A2+a_ty(i+i)A2)); 

PHI_t(i+l)  =  atan2(a_ty(i+l),a„tx(i+l)); 


125 


%  Calculate  the  components  of  the  target’s  acceleration  normal  to  the  LOS 
a_t_pitch(i+l)  =  -a_t(i+l)*cos(LAMBDA(2,i+l)-PHI_t(i+l)) 
*sin(LAMBDA(l,i+l)-TETA_t(i+l)); 

a_t_yaw(i+l)  =  -a_t(i+l)*cosOETA_t(i+l))*sin(LAMBDA(2,i+l)-PHI_t(i+l)); 
%  Update  guidance  system  input 
UgsO.i+l)  =  N*(Vc*[Xsk(2,i+l);Xsk(4,i-fl)] 

+.5*[a_t_pitch(i+l);a_t_yaw(i+l)]); 

%  Update  Time/time  to  go 
TIME(i-?-l)  =  TEME(i)+dt; 

Tgo(i+l)  =  R(4,i+ 1  )/V  c_pitch(i+ 1 ); 

%  Check  for  closest  point 
if  (R(4,i)  <  R(4,i+1)), break, end 
end; 

%  Save  information  for  plotting  and  evaluation 
Rl(4*ti+1)  a  R(4,i);  %  miss  distance 
Ti(4*ti+1)  =  ti/2;%  starting  time  of  evasive  maneuver  (EM) 
tgo(4*ti+l)  =  i*dt-ti/2;  %  time  to  go  until  end  of  flight 
if  ti  ==  12.0  %  Record  information  for  a  target  that 

%  initialized  the  evasive  maneuver  6  sec 
%  after  the  missile  entered  into  the  terminal 
TGO  =  tgo(49);  %  phase  of  flight 
Xseeker  =  Xsk(:,l:i); 

Xgsys  =  Xgs(:5l:i); 
lambda_m  =  LAMBDA_m(:,l:i); 
lambda_t  =  LAMBDA_t(:,l:i); 
lambda  =  LAMBDA(:,l:i); 
gamma_m  =  GAMMA_m(:,i:i); 
gamma_t  =  GAMMA_t(:,l:i); 
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r  =  R(:,l:i); 
vm  =  Vm(l:i); 
vt  =  Vt(l:i); 

vm_pitch  =  Vm_pitch(l:i); 
vt_pitch  =  Vt_pitch(l:i); 
vm_yaw  =  Vm_yaw(l:i); 
vt_yaw  =  Vt_yaw(l:i); 
vc_pitch  =  Vc_pitch(l:i); 
vc_yaw  =  Vc_yaw(l:i); 
tGO  =  Tgo(l:i); 
a_M  =  am(:,l:i); 
a_T  =  at(:,l:i); 

A_t  =  a_t(l:i); 

A_m  =  a_m(l:i); 
time  =  TIME(l:i); 
end 

clear  R; 

R(:,l)  =  [Xt(l,l)-Xm(l,l); 

Xt(3,l)-Xm(3,l); 

Xt(5,l)-Xm(5,l); 

sqrt((Xt(l,l)-Xm(l,l))A2+(Xt(3,l)-Xm(3,l))A2+(Xt(5,l)-Xm(5,l))A2)]; 

end; 


save  ihesis2a343  R1  tgo  Ti  tGO  missile  target  TGO  Xseeker  Xgsys  lambda. 


lambda_t  lambda  gamma_m  gamma_t  r  vm  vt  vm_pitch  vm_yaw  vt_pitch 
vt_yaw  vc_pitch  vc_yaw  a„M  a_T  time  A_t  A_m 
PLOTS 


%  Miss  distance  information 

plot(Ti,Rl),titie(‘MISS  DISTANCE  vs  INITIAL  TIME,  APROPNAV’) 


m 
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xiabel(‘ INITIAL  TIME  -  SEC’),ylabel(‘MISS  -  FEET’) 

print  -dps  Rlaal 

Ipstoepsi  Rlaal.ps  Rlaal. epsi 

pause.clg 

plot(tgo,Rl),title(‘MISS  DISTANCE  vs  TIME  TO  GO,  APROPNAV’) 

xlabel(‘TIME  TO  GO  -  SEC’).ylabel(‘MISS  -  FEET’) 

print  -dps  Rlbal 

Ipstoepsi  Rlbal.ps  Rlbal. epsi 

pause.clg 

%  Missile  acceleration  information 


plot(tirne,A_m),title(‘MISSILE  ACCELERATION  MAGNITUDE  vs  TIME, 
APROPNAV’) 

xlabel(‘TIME  -  SEC’),ylabel(‘FEET/SECA2’) 
print  -dps  A_mal 
Ipstoepsi  A_mal.ps  A_mal. epsi 
pause.clg 

plot(time,Xgsys(l,:)),title(‘MISSILE  PITCH  ACCELERATION  vs  TIME, 

APROPNAV’) 

xlabel(‘TlME  -  SEC’),ylabel(‘FEET/SECA2’) 


print  -dps  Xgsyslal 

Ipstoepsi  Xgsyslal.ps  Xgsyslal. epsi 

pause.clg 


plot(time,Xgsys(2,:)),titlc(‘MISSILE  YAW  ACCELERATION  vs  TIME, 


APROPNAV’) 

xlabel(‘TIME  -  SEC’).ylabel(‘FEET/SECA2’) 


print  -dps  Xgsys2al 


Ipstoepsi  Xgsys2al.ps  Xgsys2al.epsi 
pause.clg 
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%  Target  acceleration  information 

pIot(time,A_t),title(‘TARGET  ACCELERATION  MAGNITUDE  vs  TIME, 
APROPNAV’) 

xlabdCTIME  -  SEC’),ylabel(‘FEET/SECA2’) 

print  -dps  A_tal 

Ipstoepsi  A_tai.ps  A_tal.cpsi 

pause, clg 

%  Seeker  pitch  and  yaw  angles 

plot(time,Xseeker ( 1 )),title( ‘ SEEKEK  PITCH  ANGLE  vs  TIME, APROPNAV’) 

xlabel(‘TIME  -  SEC’),ylabel(‘RAD’) 

print  -dps  Xseekerlal 

Ipstoepsi  Xseekerlal.ps  Xseekerlal.epsi 

pause.clg 

plot(time,Xseeker(3,:)),title(‘SEEKER  YAW  ANGLE  vs  TIME,  APROPNAV’) 

xlabel(‘TIME  -  SEC’),ylabel(‘RAD’) 

print  -dps  Xsceker2al 

Ipstoepsi  Xseeker2al.ps  Xseeker2al.epsi 

pause, clg 

plot(time,Xseeker(2,:)),title(‘SEEKER  PITCH  ANGLE  RATE  vs  TIME. 

APROPNAV’) 

xlabel(‘TIME  -  SEC’),ylabel(‘RAD/SEC’) 
print  -dps  Xseeker3al 

lycawpai  /vaccine*  Jdi.epSi 

pause.clg 

plot(time,Xseeker(4,:)),title(‘SEEKER  YAW  ANGLE  RATE  vs  TIME, 

APROPNAV’) 

xiabelCTIME  -  SEC’),ylabel(‘RAD/SEC’) 
print  -dps  Xseeker4al 
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ipstoepsi  Xseeker4al.ps  Xseeker4al.epsi 
pause, cig 

%  Range  information 

plot(time,r(4,:)),titie(‘RANGE  vs  TIME,  APROPNAV’) 

xlabel(‘TJME  -  SEC),ylabel(*FEET’) 

print  -dps  ral 

ipstoepsi  ral.ps  ral.easi 

pause, clg 

%  Missile  velocity  information 

plot(time,vm),title(‘MISSILE  SPEED  vs  TIME,  APROPNAV’) 

xlabeK'TTME  -  SEC'),ylabel(‘FEET/SEC’) 

print  -dps  vrnal 

ipstoepsi  vmal.ps  vmal.epsi 

pause.clg 

%  Target  velocity  information 

plot(time,vt),titlc( ‘TARGET  SPEED  vs  TIME,  APROPNAV’) 

xlabel(‘TTME  -  SEC’),ylabel(‘ FEET/S  EC’) 

print  -dps  vtal 

ipstoepsi  vtal.ps  vtal.epsi 

pause.clg 

%  Closing  velocity  information 

plot(time,vc„pitch),title( 1  PITCH  CLOSING  SPEED,APROFNAV’) 

xlabel(‘TIME  -  SEC’),ylabei(‘FEET/SEC’) 

print  -dps  vclal 

ipstoepsi  vclal.ps  vclal. epsi 

pause.clg 

plot(time,vc_yaw),title(‘YAW  CLOSING  SPEED  vs  TIME,  APROPNAV’) 
xlabel(‘TIME  -  SEC’),ylatel(‘FEET/SEC’) 
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print  -dps  vc2al 
(pstoepsi  vc2al.ps  vc2al.epsi 
pause, clg 
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APPENDIX  C  -  MISSILE/TARGET  THREE  DIMENSIONAL 
SIMULATION  USING  OPTIMAL  GUIDANCE 


%  Written  by:  Rui  Manuel  Alves  Francisco 
°!c  Date:  09  December  1992 

%  This  Program  simulates  the  terminal  phase  of  a  3-D  missile/target 
%  engagement  using  optimal  guidance. 

clear 

clg 

%  DEFINE  CONSTANTS 
dt  =  .01;  %  Sampling  time 
Tf  =  100;  %  maximum  simulation  time 
kmax  =  Tf/dt+1; 

%  DEFINE  STATE  EQUATIONS 
%  Missile 

%  Xm  =  [xm  =  missile's  x  coordinate 
%  xmd  =  missile's  speed  (x  coordinate) 

%  ym  =  missile's  y  coordinate 

%  ymd  =  missile's  speed  (y  coordinate) 

%  zm  =  missile's  z  coordinate 

%  zmd  =  missile's  speed  (z  coordinate)] 

Am  =  [0  1  0  0  0  0 
000000 
000100 
000000 
000001 
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000000]; 
Bm  =  [000 
100 
000 
0  10 
000 
0  0  1]; 

%  Target 


% 

Xt  = [xt  = 

target’s  x  coordinate 

% 

xtd  = 

=  target’s  speed  (x  coordinate) 

% 

yt  = 

:  target’s  y  coordinate 

% 

ytd 

=  target’s  speed  (y  coordinate) 

% 

zt  = 

target’s  z  coordinate 

% 

ztd 

=  target’s  speed  (z  coordinate)] 

At 

=  [0  100 

00 

0000 

00 

000  1 

00 

0000 

00 

0000 

01 

0000 

00]; 

Bt 

=  [000 

100 

000 

010 

000 

0  0  1]; 

% 

Seeker  (Radar) 
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%  Xsk  =  [beta_pitch  =  seeker’s  pitch  angle 
%  betad_pitch  =  seeker’s  pitch  angle  rate 

%  beta_yaw  =  seeker’s  yaw  angle 

%  betad_..yaw  =  seeker’s  yaw  angle  rate] 

Ask  =  [  0  1  0  0 

-100  -20  0  0 

0  0  0  1 
0  0  -100-20]; 

Bsk  =  [  0  0 

100  0 
0  0 
0  100]; 

%  Guidance  System 

%  Xgs  =  [a_m_pitch  =  missile’s  pitch  acceleration 
%  a_m_yaw  =  missile’s  yaw  acceleration] 

Ags  =  [-10 
0-1]; 

Bgs  =  [10 

0  1]; 


%  INITIALIZE  STATE  VARIABLES  (when  the  missile  enters  into  the  terminal 


phase  of  flight) 

%  Missile 

Xm(:,l)  =  [  0 

%  The  missile  is  in  a  collision  triangle 

2828 

%  with  the  target  when  the  missile  enters  into 

0 

%  the  terminal  phase  of  flight 

1000 

0 
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47.1339]; 


%  Target 
Xt(:.l)  =  [30000 
0 
0 

1000 

500 

0]; 

%  DISCRETE  REPRESENTATION 
[PHIm.DELm]  =  c2d(Am,Bm,dt); 

[FHItJDELt]  =  c2d(At,Bt,di); 

[PHIsk.DELsk]  =  c2d(Ask,Bsk,dt); 

[PHIgs.DELgs]  =  c2d(Ags,Bgs,dt); 

%  UNE  OF  SIGHT  (LOS)  INFORMATION.  INITIAL  CONDITIONS. 
%  Missile 

%  LAMBDAjm  =  Missile’s  LOS  from  the  global  coordinate  system 
%  LAMBDAjm  =  [LAMBDA_m_pitch  =  Missile’s  pitch  LOS  angle 
%  LAMB D A_m_y aw  =  Missile’s  yaw  LOS  angle] 

LAMBDA_m(..l)  =  (atan2(Xm(5,l),sqrt(Xm(U)^+Xm(3,l)A2)); 
atan2(Xm(3,l),Xm(l,l))]; 

%  Target 

%  LAMBDAjt  =  Target’s  LOS  from  the  global  coordinate  system 
%  LAMBDA_t  =  [LAMBDA_t__pitch  =  Target’s  pitch  LOS  angle 
%  LAMBDA_t_yaw  =  Target’s  yaw  LOS  angle] 

LAMBDA_t(:,l)  =  [atan2(Xt(5,l),sqrt(Xt(l,l)^2+Xt(3,l)A2)); 

atan2(Xt(3,  l),Xt(  1,1))]; 

%  LOS  from  Missile  to  Target 
%  LAMBDA  =  LOS  from  Missile  to  Target. 
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%  LAMBDA  =  [LAMBDA_pitch  =  LOS  angle  in  pitch 
%  LAMBDA_yaw  =  LOS  angle  in  yaw]; 

LAMBDAO.l)  =  [atan2((Xt(5, 1  )-Xm(5,  l)),sqrt((Xt(  1 , 1  )-Xm(  1 ,  i  ))A2 
+(Xt(3,l)-Xm(3,l))A2)); 

atan2((Xt(3,l  )-Xm(3,  l)),abs(Xt(  1, 1)-Xm(  1 ,1)))] ; 

%  MISSILE  and  TARGET  FLIGHT  PATH  ANGLES  INFORMATION 
%  Missile 

%GAMMA_m  -  [GAMMA_m_pitch  =  Missile V,  flight  path  angle  in  pitch 
%  GAMMA_m_yaw  =  Missile’s  flight  path  angle  in  yaw] 

GAMMA_m(:,l)  =  tatan2(Xm(6,l),sqrt(Xm(2,l)A2+Xm(4,l)A2)); 
atan2(Xm(4,l),Xm(2,l))]; 

%  Target 

%GAMMA_t  =  [GAMMA_t.j)itch  =  Target’s  flight  path  angle  in  pitch 
%  G AMMA_t_Y AW  =  Target’s  flight  path  angle  in  yaw] 

GAMMA_t(:,l)  »  [atan2(Xt(6, 1  ),sqrt(Xt(2, 1)A2+Xt(4, 1  )A2)), 
atan2(Xt(4,l),Xt(2,l))]; 

%  RANGE  INFORMATION 
%  Missile 

%  Rm  =  Missile’s  range 

Rm(l)  =  sqrt(Xm(l,l)A2  +  Xm(3,l)A2  +  Xm(5,l)A2); 

%  Target 

%  Rt  =  Target’s  range 

Rt(l)  =  sqrt(Xt(l,l)A2  +  Xt(3,l)A2  +  Xt(5,l)A2); 

%  Missile/Target  relative  distance 
%  R  =  [Rmtx  =  Missth/Target  x  coordinate  range 
%  Rmty  =  Missile/Target  y  coordinate  range 

%  Rmtz  =  Missile/Target  z  coordinate  range 

%  Rmt  =  Mssile/Target  relative  distanoc(Rmt=sqrt(RmtxA2+RrntyA2+RmtzA2))] 
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R(:.l)  =  [Xt(U)-Xm(l,l) 

Xt(3,l)-Xm(3,l) 

Xt(5,l)-Xm(5,l) 

sqrt((Xt(l,  1  )-Xm(  1 , 1  ))A2+(Xt(3 , 1  )-Xm(3, 1  ))A2+(Xt(5, 1  )-Xm(5 , 1  ))A2)] , 

%  VELOCITY  INFORMATION 
%  Missile 

%  Vm  =  Missile’s  Speed 

Vm(l)  =  sqrt(Xm(2, 1  )A2+Xm(4, 1  )A2+Xm(6, 1  )A2); 

%  Target 

%  Vt  =  Target’s  Speed 

Vt(l)  =  sqrt(Xt(2, 1 )  A2+Xt(4, 1  )A2+Xt(6, 1 )  A2); 

%  Speed  along  the  pitch  and  yaw  LOS.  Pitch  and  yaw  closing  speeds 
VLpitch(l)  =  Vt(l)*cos(LAMBDA(2,l)-GAMMA_t(2,l)); 

Vm_pitch(l)  =  Vm(l)*cos(LAMBDA(2,l)-GAMMA_m(2,l)); 

Vc_pitch(l)  =  V  ra.  jjitch(  1 )  *cos(G  AMMA_.m  (1,1)  -LAMB  D  A(  1,1)) 
-Vt_pitch(l)*cos(GAMMA_t(  l ,  1  )-LAMBDA(  1.1)); 
Vt_yaw(l)  =  Vt(l)*cos(GAMMA_t(l,l)); 

Vm_yaw(l)  =  Vm(l)*cos(GAMMA_m(l,l)); 

Vc_yaw(l)  =  Vm_yaw(l)*cos(GAMMA_m(2,l)-LAMBDA(2,l)) 

-  Vt_y  a  w(  1 )  *cos(G  AMMA_t(2 , 1 )  -LAMB  D  A(2 , 1 )); 

Vc  =  [Vc_pitch(l)  0 

0  Vc_yaw(l)]; 

%  Tgo  =  Time  to  go 
Tgo(l)  =  R(4,l)/Vc_pitch(l); 

%  Optimal  guidance  coefficients 
k=Tgo(l); 

n  =  (6*kA2*(exp(-k)-l+k))/(2’,,kA3+3+6,,‘k-6*kA2-12*k,t,exp(-k)-3*exp(-2*k)); 
N  =  [n  0 
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0  n]; 

%  SEEKER  and  GUIDANCE  SYSTEM  INITIAL  CONDITIONS  and  INPUTS. 
%  Seeker 
Xsk(:.l)  =  10 
0 
0 

0]; 

Usk(:,l)  =  LAMBDA(:,1);  %  Seeker  input 
%  Guidance  System 
XgsO.l)  =  [0 
0]; 

Ugs(:J )  =  N*Vc*[Xsk(2,l)  %  Guidance  system  input 
Xsk(4,l)]; 

%  Initial  conditions  of  the  target’s  acceleration 
a_tx(l)  =0; 

a_ty(l)  =0; 

a_tz(l)  =0; 

a_t(l)  =0; 

a_t_pitch(l)  =0; 
a_t_yaw(l)  =0; 

TETA_t(  1)  =  0;  %  angle  between  the  acceleration  vector  and  the  yaw  plane 

PHI_t(  1)  =  0;  %  yaw  angle  of  the  target’s  acceleration  vector 

%  Initial  conditions  of  the  missile’s  acceleration 
a_  rn_pitch(l)  =0; 

a_m_yaw(l)  =0; 

%  TIME 

%  Time  =  Time  vector 
TTME(l)  =  0; 
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%  SIMULATE  THE  SYSTEM 
for  ti  =  0:.  25:2 1.25 
for  i  =  l:kmax-l 

%  Calculate  components  of  the  missile’s  pitch  acceleration  vector 
a_mx_pitch(i)  =  -(Xgs(l,i)*sin(LAMBDA(l,i))*cos(LAMBDA(2,i))); 
a_my_pitch(i)  =  ■(Xgs(i,i),iisin(LAMBDA(l,i))>l‘sin(LAMBDA(2,i))); 
a_mz_pitch(i)  =  Xgs(l,i)*eos(LAMBDA(l,i)); 

%  Calculate  missile’s  yaw  acceleration  vector  components 
a_mx_yaw(i)  =  -Xgs(2,i)*sin(LAMBDA(2,i)); 
a_my_yaw(i)  =  Xgs(2,i)*cos(LAMBDA(2,i)); 

%  Compute  overall  missile  acceleration 
a_mx(i)  =  a_mx_pitch(i)+a_mx  _yaw(i); 
a_my(i)  =  a.  _my_pitch(i)+a_my  _jaw(i); 
a_mz(i)  =  a_mz_pitch(i); 
am(:,<)  =  [a_mx(i) 
a_my(i) 
a_mz(i)]; 

%  target  acceleration  vector 
at(:,i)  =  [a_tx(i) 
a_ty(i) 
a_tz(i)]; 

%  Compute  magnitude  of  the  missile’s  acceleration 

A  tVI  —  mi'f/  n.  I  t  /-»  \  AO\> 

“  04*  v  -r  ^  t 

%  Generate  target’s  evasive  maneuver  (we  assume  that  these  accelerations,  along 
%  the  three  cartesian  axis,  are  estimated  using  the  missile’s  image  processing 
%  capabilities) 

if  TIME(i)  >=  ti/2  %  target  starts  evasive  maneuver 
a_tx(i+ 1)  =  3*32.2*sin(GAMMA_t(2,i)); 
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a_ty(i+l)  =  4*32.2*cos(GAMMA_t(2,i)); 
a_tz(i+l)  =  3*32.2*co$(GAMMA_t(l,i)); 
else 

a_tx(i+l)  =  0.0; 
a_ty(i+l)  =  0.0; 
a_tz(i+l)  =  0.0; 
end 

%  Compute  magnitude  of  the  target’s  acceleration 
a_t(i+l)  =  sqrt(a_tx(i+l)A2  +  a_ty(i+l)A2  +  a_tz(i-f  1)A2); 

%  Update  missile  states 

Xm(:,i+1)  =  PHIm*Xm(:,i)+DELm*am(:,i); 

%  Update  target  states 

XtC.i+1)  =  PHIt*Xt(:,i)+DELt*at(:,i); 

%  Update  flight  path  angles 

GAMMA_m(:,i+l)  =  [atan2(Xm(6,i+l),sqrt(Xm(2,i+l)A2+Xm(4.i+l)A2)); 

atan2(Xm(4,i+.  i  ),Xm(2,i+ 1 ))] ; 

GAMMA_t(:,i+l)  =  [atan2(Xt(6,i+l),sqrt(Xt(2,i+l)A2+Xt(4,i+l)A2)); 

atan2(Xt(4,i+ 1  ),Xt(2,i+ 1 ))] ; 

%  Update  seeker  states 

Xsk(:,i+1)  =  PHIsk*Xsk(:,i)+DEI..sk*Usk(:,i); 

%  Update  Guidance  System  states 
Xgs(:,i+1)  =  FHIgs*Xgs(:,i)+DELgs*Ugs(:,i); 

%  Limit  yaw  and  pitch  accelerations  to  25  g’s 
if  abs(Xgs(l,i+l))>  805.0 
Xg$(l,i+1)  =  805.0  *sign(Xgs(l,i+l)); 
end 

if  abs(Xgs(2,i+l))  >  805.0 
Xgs(2,i+1)  =  805,0  *sign(Xgs(2,i+l)); 
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end 


%  Update  LOS  angles 

LAMBDA(:,i+l)  =  [atan2((Xt(5,i+l)-Xm(5,i+l)),sqrt((Xm(l,i+l)-Xt(l,i+l))A2 
+(Xm(3,i+l)-Xt(3,i+l))A2)); 

atan2((Xt(3,i+l)-Xm(3,i+l)),(abs(Xm(l,i+l)-Xt(l,i+l))))]; 
Usk(:,i+1)  =  LAMBDA(:,i+l); 

LAMBDA_m(:,i+l)  =  [atan2(Xm(5,i+l),sqrt(Xm(l,i+l)A2+Xm(3,i+l)A2)); 

atan2(Xm(3,i+l),Xm(  1  ,i+ 1 ))] ; 

LAMB  DA _t( :  ,i+ 1 )  =  [atan2(Xt(5  ,i+  i )  ,sqrt(Xt(  1  ,i+ 1 )  A2+Xt(3  ,i+ 1 )  A2)) ; 

atan2(Xt(3,i+l),Xt(l,i+l))]; 

%  Update  Range  Information 

Rm(i+1)  =  sqrt(Xm(l,i+l)A2  +  Xm(3,i+1)A2  +  Xm(5,i+1)A2); 

Rt(i+1)  =  sqit(Xt(l,i+l)A2  +  Xt(3,i+i)A2  +  Xt(5,i+1)A2); 

R(:,i+1)  =  [Xt(l,i+1)-Xm(l,i+1); 

Xt(3,i+1)-Xm(3,i+1); 

Xt(5  ,i+ 1  )-Xm(5  ,i+ 1 ); 

sqrt((Xt(  1  ,i+ 1  )-Xm(  1  ,i+ 1  ))A2+(Xt(3,i+ 1  )-Xm(3,i+ 1 ))  A2 
+(Xt(5,i+l)-Xm(5,i+l))A2)]; 

%  Update  Velocity  Information 

Vm(i+1)  =  sqrt(Xm(2,i+l)A2+Xm(4,i+l)A2+Xm(6,i+l)A2); 

Vt(i+1)  -  sqit(Xt(2,i+l)A2+Xt(4,i+l)A2+Xt(6,i+l)A2); 

Vt_pitch(i+i)  =  Vt(i+l)*cos(LAMBDA(2,i+l)-GAMMA_t(2,i+l)); 
Vm_pitch(i+ 1 )  =  Vm(i+  l)*cos(LAMBD  A(2,i+ 1  )-GAMMA_m(2,i+ 1 )); 
Vc_pitch(i+1)  =  V  m_pitch(i+ 1 )  *cos(G  AMMA_m(  1  ,i+ 1 )-  L  AMB  D  A(  1 ,  i+ 1 )) 

-  Vt_pitch(i+ 1 )  *cos  (G  AMMA_t(  1  ,i+l)-LAMBDA(  1  ,i+ 1 )); 
Vt_yaw(i+1)  =  Vt(i+i)*cos(GAMMA_t(l,i+l)); 

Vm_yaw(i+1)  =  Vm(i+l)*cos(GAMMA_rn(l,i+l)); 

Vc_yaw(i+1)  =  Vm_yaw(i+l):i'cos(GAMMA_m(2,i+l)-LAMBDA(2,i+l)) 
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-Vt_yaw(i+l)*cos(GAMMA_t(2,i+i)-LAMBDA(2,i+l)); 

Vc  =  [Vc_pitch(i+1)  0 

0  Vc_yaw(i+1)]; 

%  Time  to  go 

Tgo(i+l)  =  R(4,i+1)/Vc_pitch(i+1); 

%  Calculate  angles  of  the  target’s  acceleration 
TETA_t(i+l)  =  atan2(a_tz(i+l),sqrt(a_tx(i+l)A2+a_ty(i+l)A2)); 

PHI_t(i+l)  =  atan2(a_ty(i+ 1  ),a_tx(i+ 1 )); 

%  Calculate  the  components  of  the  target’s  acceleration  normal  to  the  LOS 
a_t_pitch(i+l)  =  -a_tCi+l)*cos(LAMBDA(21i+l)-PHLt(i+l)) 

*sin(LAMBDA(  1  ,i+l)-TETA_t(i+ 1 )); 

a _t_yaw(i+ 1 )  =  •  a_t(i+ 1  )*co$(TETA_t(i+ 1  ))*sin(LAMB  DA(2,i+ 1  )-PHI_t(i+ 1 )) ; 
%  Update  optimal  guidance  coefficients 
k  =  Tgo(i+l); 

n  =  (65itkA25|'(exp(-k)-l+k))/(2*kA3-4-3+6*k-6!|ckA2-12*k*exp(-k)-3*exp('2*k)); 

N  =  [n  0 
0  n]; 

%  Components  of  the  Missile’s  acceleration  normal  to  the  pitch  and  yaw  LOS 
a_m_pitch(i+l)  =  Xgs(l,i+1); 
a_m_yaw(i+l)  =  Xgs(2,H-l); 

%  Update  guidance  system  input 
Ugs(:,i+1)  =  N*(Vc*[Xsk(2,i+l);Xsk(4,i+l)] 

4-.5*[a_t_pitch(i+l);a_t_yaw(i+l)]) 

-( l/kA2)*(exp(-k)+k- 1  )*[a_m_pitch(i+  l);a_m_yaw(i+ 1)] ); 

%  Update  Time 
TIME(i+l)  =TIME(i)+dt; 

%  Check  for  closest  point 
if  (R(4,i)  <  R(4,i4-l)),break,end 


end; 

%  Save  information  for  plotting  and  evaluation 
Rl(4*ti+1)  =  R(4,i);  %  miss  distance 
Ti(4*ti+1)  =  ti/2;%  starting  time  of  evasive  maneuver  (EM) 
tgo(4*ti+l)  =  i*dt-ti/2;  %  time  to  go  until  end  of  flight 
if  ti  ==  12.0  %  Record  information  for  a  target  that 

%  initialized  the  evasive  maneuver  6  sec  after 
%  the  missile  entered  into  the  terminal  phase 
TGO  =  tgo(49);  %  of  flight 
Xseeker  =  Xsk(:,l:i); 

Xgsys  =  Xgs(:,l:i); 
lambda_m  =  LAMBDA_m(:,l:i); 
lambda_t  =  LAMBDA_t(:,l:i); 
lambda  =  LAMB  DA(:,l:i); 
gamma_m  =  GAMMA_m(:,l:i); 
gamma_t  =  GAMMA_t(:,l:i); 
r  =  R(:,l:i); 
vm  =  Vm(l;i); 
vt  =  Vt(l:i); 

vm_pitch  =  Vm_pitch(l;i); 
vt_pitch  =  Vt_pitch(l:i); 
vm_yaw  =  Vm._yaw(l:i); 
vt_yaw  =  Vt_yaw(l:i); 
vc_pitch  =  Vc_pitch(l:i); 
vc_yaw  =  Vc_yaw(l:i); 
tGO  =  Tgo(i:i); 
a_M  -•  am(:,l:i); 
a_T  =  at(:,l:i); 
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A_t  =  a_t(l:i); 

A_m  =  a_m(l:i); 

time  =  TIME(l:i);  ^ 

end 

clear  R;  * 

R(:,i)  =  [Xt(l,l)-Xm(l,l); 

Xt(3,l)-Xm(3,l); 

Xt(5,l)-Xm(5,l); 

sqrt((Xt(  1 , 1  )-Xm(  1 , 1  ))A2+(Xt(3, 1)-Xm(3, 1  ))A2+(Xt(5 ,  1)-Xm(5 , 1  ))A2)]; 
end; 

save  thesis3o343  R1  tgo  Ti  tGO  missile  target  TGO  Xseeker  Xgsys  lambda_m 
lambda_t  lambda  gamrna_m  gamma_t  r  vm  vt  vm_pitch  vm_yaw  vt_pitch 
vt_yaw  vc_pitch  vc_yaw  a_M  a_T  time  A_.t  A_m 
PLOTS 

%  Miss  distance  information 

piot(Ti,R  1  ),title(‘MISS  DISTANCE  vs. INITIAL  TIME,  OPTIMAL’) 

xlabel( ‘INITIAL  TIME  -  SEC’),ylabel(‘MSS  -  FEET’) 

print  -dps  Rlaol 

Ipstoepsi  Rlaol.ps  Rlaol. epsi 

pause,  clg 

plot(tgo,Rl),title(‘MISS  DISTANCE  vs  TIME  TO  GO,  OPTIMAL’) 

xlabel(‘TIME  TO  GO  -  SEC’),ylabel(‘MISS  -  FEET’) 

print  -dps  Rlbol 

Ipstoepsi  Rlbol.ps  Rlbol.epsi 

pause, clg 

%  Missile  acceleration  information 

plot( time, A_m) ,title( ‘ MISS ILE  ACCELERATION  MAGNITUDE  vS  TIME, 

OPTIMAL’) 
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xlabel(‘TIME  -  S  EC  ’ )  ,y  labe  1(  ‘  FEET/S  EC A2  ’ ) 
print  -dps  A_mol 
Ipstoepsi  A_mol.ps  A_mol.epsi 
pause, clg 

plot(time,Xgsys(l.:)).title(‘MISSILE  PITCH  ACCELERATION  vs  TIME, 

OPTIMAL’) 

xlabelOTEME  -  SEC’),ylabel(‘FEET/SECA2’) 
print  -dps  Xgsyslol 
Ipstoepsi  Xgsyslol.ps  Xgsyslol.epsi 
pause.clg 

plot(time,Xgsys(2,:)),title(‘MISSILE  YAW  ACCELERATION  vs  TIME, 

OPTIMAL') 

xlabei(‘TIME  -  S  EC '), y  labe  1(‘ FEET/S  EC  A2’) 
print  -dps  Xgsys2ol 
Ipstoepsi  Xgsys2ol.ps  Xgsys2ol,epsi 
pause.clg 

%  Target  acceleration  information 

plot(time, A_t) ,title(‘TARGET  ACCELERATION  MAGNITUDE  vs  TIME, 
OFITMAL’) 

xlabel(‘TIME  -  SEC’),ylabel(‘FEET/SECA2’) 

print  -dps  A_tol 

Ipstoepsi  A_tol.ps  A_tol.epsi 

pause.clg 

%  Seeker  pitch  and  yaw  angles 

plot(time,Xseeker(  1 ,:)), title(‘ S EEKER  PITCH  ANGLE  vs  TIME,  OPTIMAL’) 

xlabcl(‘TIME  -  SEC’),ylabel(‘RAD’) 

print  -dps  Xseekerlol 

Ipstoepsi  Xseekerlol.ps  Xseekerlol.eps: 
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pause, clg 

plot(time,Xseeker(3,:)),title(‘SEEKER  YAW  ANGLE  vs  TIME,  OPTIMAL’) 

xlabelCTIME  -  SEC’).ylabel(‘RAD’) 

print  -dps  Xseeker2ol 

Ipstoepsi  Xseeker2ol.ps  Xseeker2ol.epsi 

pause, clg 

pIot(time,Xseeker(2,:)),title(‘SEEKER  PITCH  ANGLE  RATE  vs  TIME, 

OPTIMAL’) 

xlabel(‘TIME  -  SEC’),ylabel(‘RAD/SEC’) 
print  -dps  Xseeker3ol 
Ipstoepsi  Xseeker3ol.ps  Xseeker3ol.epsi 
pause, clg 

plot(time,Xseeker(4,:)),title(‘SEEKER  YAW  ANGLE  RATE  vs  TIME, 

OPTIMAL’) 

xlabelCTIME  -  SEC’),ylabel(‘RAD/SEC’) 
print  -dps  Xseeker4ol 
Ipstoepsi  Xseeker4ol.ps  Xseeker4ol.epsi 
pause.clg 

%  Range  information 

plot(time,r(4,:)),title(‘RANGE  vs  TIME,  OPTIMAL’) 

xlabelCTIME  -  SEC’),ylabel(‘FEET’) 

print  -dps  rol 

Ipstoepsi  rol.ps  rol.easi 

pause, clg 

%  Missile  velocity  information 

plot(time,vm),title(‘MISSILE  SPEED  vs  TIME,  OPTIMAL’) 
xlabelCTIME  -  SEC’),ylabcl(‘FEET/SEC’) 
print  -dps  vmol 
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Ipstoepsi  vmol.ps  vmol.epsi 
pause, clg 

%  Target  velocity  information 

plot(time,vt),title(‘TARGET  SPEED  vs  TIME,  OPTIMAL’) 

xlabeK'TIME  -  SEC’),ylabel(‘FEET/SEC’) 

print  -dps  vtol 

!pstoepci  vtol.ps  vtol.epsi 

pause, clg 

%  Closing  velocity  information 

plot(time.vc_pitch),title(‘PITCH  CLOSING  SPEED,  OPTIMAL’) 

xlabeK’TIME  -  SEC’),ylabel(‘FEET/SEC’) 

print  -dps  vclol 

Ipstoepsi  vclol.ps  vclol. epsi 

pause.clg 

plot( time, vc_yaw),title(‘ YAW  CLOSING  SPEED  vs  TIME,  OPTIMAL’) 

xiabel(‘TIME  -  SEC’),yiabel('FEET/SEC’) 

print  -dps  vc2ol 

Ipstoepsi  vc2ol.ps  vc2ol.epsi 

pause.clg 
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